function myODE()
evalin('base','clear')
clc
close all
tspan=inputdlg('Enter timespan:','tspan',1,{'enter time in seconds'});
tspan=str2double(tspan{1});
t = linspace(0,tspan, 100*tspan);
x0 = [0; 0; eps; eps]; % x0, y0, vx0, vy0
% equation blows up with zero initial velocities
eqn1=inputdlg('Enter in xddot equation:','xddot',1,{'x double dot ='});
eqn1=eqn1{1};
eqn1=['@(t,xdot,ydot,a,g,uk,m)',eqn1];
fun1=str2func(eqn1);
eqn2=inputdlg('Enter in yddot equation:','yddot',1,{'y double dot ='});
eqn2=eqn2{1};
eqn2=['@(t,xdot,ydot,b,g,uk,m)',eqn2];
fun2=str2func(eqn2);
opts = odeset('Reltol',1e-13,'AbsTol',1e-14,'Stats','on');
fun=@(t,x)xddot(t,x,fun1,fun2);
[t, x] = ode113(fun, t, x0, opts);
figure;
subplot(1,2,1); plot(x(:,3),x(:,4));
xlabel('v_x', 'fontsize', 20); ylabel('v_y', 'fontsize', 20);
title('Velocity in the plane parameterized by t')
subplot(1,2,2); plot(x(:,1),x(:,2));
xlabel('x', 'fontsize', 20); ylabel('y', 'fontsize', 20);
title('Position in the plane parameterized by t')
end
function out = xddot(t,x,fun1,fun2)
m = 1;
g = 9.81;
uk = .3;
a = 10;
b = 10;
xdot=x(3);
ydot=x(4);
out = [0; 0; 0; 0];
out(1) = xdot;
out(2) = ydot;
out(3) = fun1(t,xdot,ydot,a,g,uk,m);
out(4) = fun2(t,xdot,ydot,b,g,uk,m);
end