Using Implicit Euler Method with Newton-Raphson method

So I'm following this algorithm to write a code on implicit euler method

and here is my attempt

function y = imp_euler(f,f_y,t0,T,y0,h,tol,N)
t = t0:h:T;
n = length(t);
y = zeros(n,1);
y(1) = y0;
for k = 1:n-1
    g = @(z) z - y(k) - h*f(t(k+1),z);
    gp = @(z) 1 - h*f_y(t(k+1),z) ;
    y(k+1) = newton(f,f_y,y(k),tol,N);
end
end

where

function sol=newton(f,fp,x0,tol,N)
    i=0;
    sol = zeros(N,2);
    fc=abs(f(x0));
    while (fc>tol)
        xc = x0 - (f(x0)/fp(x0));
        fc=abs(f(xc));
        x0 = xc;
        i=i+1;
        sol(i,:) = [i; x0];
        if (i>N)
            fprintf('Method failed after %d iterations. \n',N);
            break
        end
    end
    sol = sol(any(sol,2),:);
end

Unfortunately, my code does not work for some reason. Could anybody guide me on how to code this? Comments are appreciated.

[Torsten]

You forgot to include the call to imp_euler.

[Eugene Miller]

could you please elaborate?

[Torsten]

I mean: The code you provided doesn't work without the knowledge of f, fy, t0,T,y0,h,tol and N (I think n).

[Eugene Miller]

its a function

[Torsten]

How did you come the conclusion that the code does not work if you didn't test it ? Because for testing, inputs of the mentionned variables had to be provided.

[Eugene Miller]

i tested it with this

f = @(t, y) -20*t*y^2;

f_y = @(t, y) -40*t*y;

t0 = 0;

T = 1;

y0 = 1;

h = 0.2;

tol = 1e-8;

N = 100;

[Eugene Miller]

and the answer is

y =

0.2

(but it should be a vector)

Answer 1

[Torsten]

f = @(t,y) -20*t*y^2;
f_y = @(t,y) -40*t*y;
t0 = 0;
T = 1;
y0 = 1;
h = 0.01;
tol = 1e-8;
N = 100;
y = imp_euler(f,f_y,t0,T,y0,h,tol,N)
plot(t0:h:T,y)
hold on
plot(t0:h:T,1./(1+10*(t0:h:T).^2))


function y = imp_euler(f,f_y,t0,T,y0,h,tol,N)
t = t0:h:T;
n = length(t);
y = zeros(n,1);
y(1) = y0;
for k = 1:n-1
    g = @(z) z - y(k) - h*f(t(k+1),z);
    gp = @(z) 1 - h*f_y(t(k+1),z) ;
    disp('Vor Newton call.')
    y(k+1) = newton(g,gp,y(k),tol,N);
    disp('nach Newton Call.')
end
end

function sol=newton(f,fp,x0,tol,N)
i=0;
sol = zeros(N,2);
fc=abs(f(x0));
while fc > tol
    xc = x0 - (f(x0)/fp(x0));
    fc=abs(f(xc));
    x0 = xc;
    i=i+1;
    if (i>N)
        fprintf('Method failed after %d iterations. \n',N);
        break
    end
end
sol = x0;
end

[Haseeb Hashim]

@Torsten Man thanks. You are a legend. I was looking for this code for so long