I want to solve the following Cauchy problem \begin{equation} y' = y^2 + \frac{t^4 - 6t^3 + 12t^2 - 14t + 9}{(1+t)^2} \end{equation}
with initial condition: $y(0) = 2$ for $t \in [0,1.6]$ using a 3 steps Adam-Moulton method implemented in Matlab as follows:
function [t,u] = AM3(fun,t0,T,y0,N)
h = (T-t0)/N; % integration step
t = t0:h:T; %time mesh
[t(1:3),u(1:3)] = RK4(fun,t0,t0+2*h,y0,2); % first 3 steps
for n=1:N-2
fn = feval(fun,t(n),u(n));
fn1 = feval(fun,t(n+1),u(n+1));
fn2 = feval(fun,t(n+2),u(n+2));
F = @(z) z - u(n+2) - (h/24)*(9*feval(fun,t(n+3),z) + 19*fn2 - 5*fn1 + fn);
z0 = u(n+1);
z1 = u(n+2);
dF = @(z) 1 - (h/24)*9*(2*z);
u(n+3) = newton(F,dF,z0);
end
endfunction
N is the number of integration nodes used (i.e. the size of the time mesh). The function returns a discretized solution of the problem in the array $u$. Since the method is implicit at every iteration I use Newton's method written as
function y = newton(f,df,x)
fx = f(x);
tol = 1e-10;
itermax = 1e3;
iter = 0;
while abs(fx) > tol && iter < itermax
x = x - f(x)/df(x);
fx = f(x);
iter++;
endwhile
y=x;
My problem is that it seems that the equation is very ill-conditioned since even with a tiny integration step $h$ what I obtain is something like this:
the problem persisting with much bigger $N$.
Despite the equation being pretty ugly it has en exact solution: $y(t) = \frac{(t-1)(t-2)}{t+1}$ which is relatively simple.
Can you help me please getting around this? I double checked the algorithms and they seems fine; I also tried to change Newton's method with some other (secant, fixed point, etc) but no novelty occured.
I also tried changing altoghether the method used with an Adam-Bashforth or Runge-Kutta one but to no avail.