Midpoint integration with for loop

Hi,

Computing a numerical integration with the Midpoint Method I'm struggling with the output of my function, while with Euler I got the expected result, with this method I'm getting a different graphical solution far away than the analytical solution (and the euler's one). What could be wrongly specified?

Thanks for your help.

% Midpoint method
% General information
L = 0.61;  % [m]
G = 9.81;  % gamma [m/s^2]

% Sampling interval
h = 0.1;  
%h = 0.05;
%h = 0.01;

ti = 0; % Initial time
tf = 4; % Final time
t  = ti:h:tf;  % Time vector

function [md] = midpoint(t,L,G,h)
% Initial Value problem
md = zeros(2, length(t));
md(1,1) = pi/40;

for i = 1:length(t)-1
    pmd = md(2,i) + (h/2) + (-G/L)*sin(md(1,i)) + (h/2) * [md(2,i); (-G/L)*sin(md(1,i))];
    md(:,i+1) = md(:,i) + h * pmd;
end
end

[Hugo Hernández Hernández]

I have updated my code, got a better result but still the same problem, this comes with "pmd" where as I include an array inside another array there comes an error, I know that this is not possible but, I don't find another solution that generates this midpoint condition, and also applies to Runge-Kutta for 4th order.