Matlab, numerical integration and two coupled nonlinear second order ode using Modified Euler, Runga-Kutta, ode45

Question

My problem is : $$ h''(T) = 1/2 \int_0^1 \left\{1 - \frac{(h(T)+(1-b)\,⋅\theta(T))^2}{(h(T)+(x-b)⋅\theta(T)-K⋅\sin(\pi⋅x))^2}\right\} dx $$ $$ \theta''(T) = \frac{1}{2⋅\Gamma} \int_0^1 \left\{(x-b)⋅\left(1 - \frac{(h(T)+(1-b)⋅\theta(T))^2}{(h(T)+(x-b)⋅\theta(T)-K⋅\sin(\pi⋅x))^2}\right)\right\} dx $$

I need to find $h$ and $\theta$ which are unknown functions of time $T$. Initial conditions : $h(0) = \kappa⋅\sin (\pi⋅b)$, $\theta(0) = \kappa⋅\pi⋅\cos(\pi⋅b)$, $h'(0)=0$, $\theta'(0)=0$.

I have tried these three methods in Matlab to compare each solution from each solver, namely modified Euler, Matlab ode45 and Runga-Kutta 4th. For integrations I have used simpsons algorithm. The plot of results for h and theta from different solvers : modified euler, ode45 and runga-kutta are very much different. Here is the respective code:

function euler
% A second order differential equation
Tsim = 1;                                            % simulate over Tsim seconds
dt = 0.01;                                           % step size
N= Tsim/dt;                                          % number of steps to simulate
x=zeros(4,N);
t = zeros(1,N);
%--------------------------------------------------------------%|
b = 0.4;                                                       %|
kappa = 1;                                                     %|
M = 1;                                                         %|
g = .01;                                                       %|
Gamma = 0.2;                                                   %|
% Shape function               % F = kappa*x*(1-x);            %|
h0 = kappa*sin(pi*b);          % h0 = F(b);                    %|
theta0 = kappa*pi*cos(pi*b);   % htheta0 = F'(b);              %|
%--------------------------------------------------------------%|

% x is the solution vector x(1,:) = h, x(2,:) = h',
% x(3,:) = theta, x(4,:) = theta'
% For h
x(1,1)= h0;        % h(t = 0) ;
x(2,1)= 0;         % dh/dt (t = 0) = h1;
% For theta
x(3,1)= theta0;    % theta(t = 0);
x(4,1)= 0;         % dtheta/dt(t = 0) = theta1;
for k=1:N-1
    % For time t
    t(k+1)= t(k) + dt;
    xnew=x(:,k)+ dt *MassBalance(t(k),x(:,k));
    x(:,k+1) = x(:,k) + dt/2* (MassBalance(t(k),x(:,k)) + MassBalance(t(k),xnew)); %Modified Euler algorithm
end
%Compare with ode45
x0 = [h0; 0; theta0; 0];
[T,Y] = ode45(@MassBalance,[0 Tsim],x0);   % where Y is the solution vector
[aa, bb]= RK4(h0, theta0, 0, 0, Tsim, N);
end

function dx= MassBalance(t,w)
% use a structure to store the parameters
%--------------------------------------------------------------%|
b = 0.4;                                                       %|
kappa = 1;                                                     %|
M = 1;                                                         %|
g = .01;                                                       %|
Gamma = 0.2;                                                   %|
% Shape function               % F = kappa*x*(1-x);            %|
h0 = kappa*sin(pi*b);          % h0 = F(b);                    %|
theta0 = kappa*pi*cos(pi*b);   % htheta0 = F'(b);              %|
%--------------------------------------------------------------%|
% Equations of motion for the body
dx = zeros(4,1);

xx = (0:0.003:1);

int_func1 = 1/2*(1 - ((w(1)+ (1-b)*w(3))./(w(1)+(xx-b).*w(3)-kappa.*sin(pi.*xx))).^2);
int_func2 = (xx-b).*(1/2*(1 - ((w(1)+ (1-b)*w(3))./(w(1)+(xx-b).*w(3)-kappa.*sin(pi.*xx))).^2));

% plot(xx, integral_func)
I1 = simpsons(int_func1,0,1,[]);
I2 = simpsons(int_func2,0,1,[]);
dx(1)=w(2);
dx(2)=I1 - M*g
dx(3)=w(4);
dx(4)= 1/Gamma.* I2;
end

function [x y u w]= RK4(a0, b0, c0, d0,  t, n)
% x = h; y = theta; u = h'; w = theta';
%--------------------------------------------------------------%|
b = 0.4;                                                       %|
kappa = 1;                                                     %|
M = 1;                                                         %|
g = .01;                                                       %|
Gamma = 0.2;                                                   %|
% Shape function               % F = kappa*x*(1-x);            %|
h0 = kappa*sin(pi*b);          % h0 = F(b);                    %|
theta0 = kappa*pi*cos(pi*b);   % htheta0 = F'(b);              %|
%--------------------------------------------------------------%|



% Runga-Kutta scheme for 2 coupled nonlinear second order system
%-----------------------------------------------------
h = t / n;
x = [ a0 ]; y = [ b0 ]; u = [ c0 ]; w = [ d0 ];

xx = (0:0.003:1);

int_func1 = 1/2*(1 - ((x+ (1-b)*y)./(x+(xx-b).*y-kappa.*sin(pi.*xx))).^2)
int_func2 = (xx-b).*(1/2*(1 - ((x+ (1-b)*y)./(x+(xx-b).*y-kappa.*sin(pi.*xx))).^2))

I1 = simpsons(int_func1,0,1,[]);
I2 = simpsons(int_func2,0,1,[]);


f1 = @(x, y, u, w) I1 - M*g;
f2 = @(x, y, u, w) I2;


for i = 1 : n
    kx1 = h * u(i);
    ky1 = h * w(i);
    ku1 = h * f1(x(i), y(i), u(i), w(i));
    kw1 = h * f2(x(i), y(i), u(i), w(i));

    kx2 = h * (u(i) + ku1/2);
    ky2 = h * (w(i) + kw1/2);
    ku2 = h * f1(x(i)+kx1/2, y(i)+ky1/2, u(i)+ku1/2, w(i)+kw1/2);
    kw2 = h * f2(x(i)+kx1/2, y(i)+ky1/2, u(i)+ku1/2, w(i)+kw1/2);

    kx3 = h * (u(i) + ku2/2);
    ky3 = h * (w(i) + kw2/2);
    ku3 = h * f1(x(i)+kx2/2, y(i)+ky2/2, u(i)+ku2/2, w(i)+kw2/2);
    kw3 = h * f2(x(i)+kx2/2, y(i)+ky2/2, u(i)+ku2/2, w(i)+kw2/2);

    kx4 = h * (u(i) + ku3);
    ky4 = h * (w(i) + kw3);
    ku4 = h * f1(x(i)+kx3, y(i)+ky3, u(i)+ku3, w(i)+kw3);
    kw4 = h * f2(x(i)+kx3, y(i)+ky3, u(i)+ku3, w(i)+kw3);

    x(i+1) = x(i) + (kx1 + 2*kx2 + 2*kx3 + kx4) / 6;
    y(i+1) = y(i) + (ky1 + 2*ky2 + 2*ky3 + ky4) / 6;
    u(i+1) = u(i) + (ku1 + 2*ku2 + 2*ku3 + ku4) / 6;
    w(i+1) = w(i) + (kw1 + 2*kw2 + 2*kw3 + kw4) / 6;
end
end

Any idea please, why I am getting very different results from these three solvers.