help fixing the Heuns method

So I need some help with this heuns method, There is something wrong in my code but I don't know where, since when i plot it the Euler method i have is more close to the function than the Heun

f: function, i.e. right side of equation y '= f (x, y)

t0 and y0: initial values such that y (x0) = y0

no_steg: number of Heun steps The Heun function should calculate

lengde: number that indicates the length of an HEun step

t0 = 0; %start value of x
y0 = 2; %start value of y
no_steg = 3;
lengde = 0.2;
t = 0:lengde:no_steg;
f= @(t,y) -2*t*y; % function 
[xverd, tilnaer] = Heun_met (f, t0, y0, no_steg, lengde)
plot(xverd, tilnaer)
hold on 
y1= y0;
for s=1: length(t)-1 
    t0(s+1) = t0(s) + lengde;
    y1(s+1)= y1(s) + lengde*f(t0(s), y1(s));    %Beregner et euler steg for startverdi problemet
end
plot(xverd, y1)
y2 = 2*exp(-t.^2);
plot(xverd, y2)
hold off

function [xverd, tilnaer] = Heun_met (f, t0, y0, no_steg, lengde)
t = 0:lengde:no_steg;
y1 = y0;

for s=1: length(t)-1
    
    t0(s+1) = t0(s) + lengde;
    y1(s+1)= y1(s) + lengde*f(t0(s), y1(s));    %Beregner et euler steg for startverdi problemet
    y0(s+1) = y0(s) + (lengde/2) * (f(t0(s), y0(s)) + (f(t0(s+1) + lengde, y1(s+1))));    %Beregner heun steg
    
end
tilnaer = y0;
xverd = t0;
end 

Answer 1

[Jim Riggs]

Heun's method is otherwise known as "explicit trapezoidal method", "improved Euler's method" or "modified Euler's method".

The way you have it,

t0(s+1) = t0(s) + lengde, 
%therefore, 
t0(s+1) + lengde = t0(s) + 2*lengde

therefore, "lengde" is being added twice, which is incorrect

[see comment, below]

[Jim Riggs]

function [t, y] = Heun_met (f, t0, y0, no_steg, lengde)
t = t0:lengde:no_steg;
y = 0.*t;
y(1) = y0;
for s=1: length(t)-1
    
    k1 = lengde * f(t(s), y(s));  
    k2 = lengde * f(t(s) + lengde, y(s) + k1)
    y(s+1) = y(s) +  (k1 + k2)/2; 
    
end
end

[Jim Riggs]

The other way to code it is:

...
for s=1: length(t)-1
    
    k1 =  f(t(s), y(s));  
    k2 =  f(t(s) + lengde, y(s) + lengde*k1)
    y(s+1) = y(s) +  lengde * (k1 + k2)/2;  
    
end
...
    

This is equivalent to the one, above.