
function msv()

global a
a = 0;

[t,y] = ode45(@f,[0,1],1);
[tmsv,ymsv] = msvint(@f,[0,1],1,25); % N =11 h=0.1

plot(t,abs(y),'r-',tmsv,abs(ymsv),'bo-'); 
xlabel('t')
ylabel('y')
legend('ode45','HoMsv')

end

function [t,y] = msvint(fh,tspan,y0,N)
global a

t0 = tspan(1);
T = tspan(2);
h = (T-t0)/N;

t = zeros(N+1,1);
y = zeros(N+1,1);

% Startwerte
t(1) = t0;
t(2) = t0+h;
y(1) = y0;
%y(2) = y0 + h*f(t(1),y(1));
y(2) = exp(a*h)*y(1)+1e-12;
for n = 0:N-2
    t(n+3) = t(n+2) + h;
    y(n+3) = msvstep(fh,y(n+2),y(n+1),t(n+1),h);
end
end

function yn2 = msvstep(fh,yn1,yn,tn,h)
  yn2 = -4*yn1 + 5*yn + h*(4*fh(tn+h,yn1) + 2*fh(tn,yn));
end

function val = f(t,y)
global a;
val = a*y;
%val = 0;
end