Taylor expansion does not work

edit

I started with Sagemath and wanted to use it in order to get the Taylor Series of a function. However, a quick check shows that the coefficients are not correct. If I use series() instead of taylor(), then the output looks okay.

Did I make something wrong? Does Sagemath have some problems as the function contains several variables?

Here is a minimal example which reproduce the problem:

reset()

var('exx eyy x kx')

p1 = kz_s = sqrt( eyy - kx^2)

p2 = sqrt(1 - kx^2)

rs = (p1 - p2) / (p1 + p2)

rs = rs(eyy = exx + x)

rs = rs.subs(eyy = exx + x)

ser = rs.taylor(x, 0, 1)

p1 = plot(rs(exx = 4, x = 0.0) , (x, 0, .2), color='red' )

p2 = plot(ser(exx = 4, x = 0.0) , (x, 0, .2) )

show(p1 + p2)