Note that there is a discrepancy between your code and your graphics for the second ode. I took the version from the graphics ( which worked better for ode15i :-) )
y0 = [0 0 0 0];
yp0 = [0 0 0 0];
t0 = 0;
tspan = [t0,0.25];
[y0_new,yp0_new] = decic(@f,t0,y0,[1 1 1 1],yp0,[0 0 0 0])
sol = ode15i(@f,tspan,y0_new,yp0_new,odeset('RelTol',1e-8,'AbsTol',1e-8));
figure(1)
plot(sol.x,[sol.y(1,:);sol.y(2,:)])
[y,yp] = deval(sol,sol.x);
for i = 1:numel(sol.x)
res(i,:) = f(sol.x(i),y(:,i),yp(:,i));
end
figure(2)
plot(sol.x,[res(:,3),res(:,4)])
function res = f(t,y,yp)
I = 0.000122144; %kgm^2
m = 0.19744; %kg
g = 9.80665; %m/s^2
R = 0.035; %m
s = 0.090757; %rad
r = -0.011515; %m
M = 0.029244; %kg
k = 15.36;
A = 0.011065331; %m^2
res = zeros(4,1);
res(1) = yp(1) - y(3);
res(2) = yp(2) - y(4);
res(3) = I*yp(3) - (-m * sqrt(g^2 + R^2 * yp(4)^2 - 2 * g * R * yp(4) * sin(s)) * r * sin(s + y(1) - atan(R * yp(4) * cos(s) / (g - R * yp(4) * sin(s)))) - k * R * A * (y(3) - y(4)));
res(4) = 2 * M * R * yp(4) -( (M + m) * g * sin(s) - m * r * (sin(y(1)) * (y(3)^2) - cos(y(1)) * yp(3)));
end
