Question: Reducing Run-time for solving a nonlinear ODE

Hi

I want to employ FDM to solve nonlinear ODE. Since large expressions in terms of s[1] are generated, the program lose its efficiency for N>6 (It requires long run-time). Please amend the program, if it is possible to reach more precision for N>10. Moreover, is there a Maple command to dsolve this ode?

Thanks alot

restart;

sy := 2400.:

Hp := 0.1e9:

P := -900.:

a := 20.:

c := 21.:

N := 6:

s[0] := P:

h := (c-a)/N:

Et := (r*(diff(sr(r), r))-sy)/Hp:

Er := (sy-r*(diff(sr(r), r)))/Hp:

ode := simplify(Er-Et-r*(diff(Et, r))*(1+(diff(Et, r))*r/(2+4*Et))):

ode:=subs(diff(sr(r), r, r) = (s[k+1]-2*s[k]+s[k-1])/h^2, diff(sr(r), r) = (s[k+1]-s[k-1])/(2*h), r = h*k+a, ode):

for k to N-1 do

s[k+1] := solve(ode, s[k+1])[1]

end do:

s[1] := fsolve(s[N] = -200., s[1] = -900 .. -100);

plots[pointplot]([seq([h*k+a, s[k]], k = 0 .. N)]);