Hi
Not sure if I understand the question correctly, but you basically already have the answer, and most of the code you need.
Assuming you want to implement your own version of Runge Kutta and not use the matlab ode, you just need to slightly change 'myfun':
Yh=[z0;zd0];
for i=2:length(time)
dY = myode(time(i),dt,Yh);
temp=Yh(:,i-1)+dY*dt;
Yh=cat(2,Yh,temp);
end
Just replace the line with % Euler integration by the Runge Kutta scheme of your choice. The most popular one is RK4, you'll find everything you need on wikipedia: Runge Kutta Methods One suggestion though: Instead of concatenating your result as Yh=cat(2,Yh,temp), you should preallocate the whole vector, e.g.:
Y = zeros(2,length(time));
Y(:,1) = [z0;zd0];
for i = 2:length(time)
...
Yh(:,i) = Yh(:,i-1) + ...
end
