Execute MATLAB/Octave Online

clear all
close all
clc

%Gaurav Sharma
%AERO 6330
%Hw 3

format long;

time = 0.3; %total time of computations
delt = 0.01; %step size

%Initial conditions
t(1) = 0.0;
x(1) = 1.0;
y(1) = 2.0;
z(1) = 3.0;
vx(1) = 2.0;
vy(1) = 3.0;
vz(1) = 4.0;

syms t vx vy vz x y z;
fx = symfun (6*vx-9*x, [t x y z vx vy vz]);
syms t vx vy vz x y z;
fy = symfun (6*vy-9*y, [t x y z vx vy vz]);
syms t vx vy vz x y z;
fz = symfun (6*vz-9*z, [t x y z vx vy vz]);

vxdot = fx(t, vx, vy, vz, x, y, z); %(t,r,v); %fx = 6*vx-9*x;
vydot = fy(t, vx, vy, vz, x, y, z); %(t,r,v); %fy = 6*vy-9*y;
vzdot = fz(t, vx, vy, vz, x, y, z); %(t,r,v); %fz = 6*vz-9*z;

%vdot = [vxdot vydot vzdot];
% v = [vx vy vz];
% r = [x y z];

index = 0.0;

for i=1:delt:time+1;

k1x = delt*fx(t,x,y,z,vx,vy,vz);
k1y = delt*fy(t,x,y,z,vx,vy,vz);
k1z = delt*fz(t,x,y,z,vx,vy,vz);

k2x = delt*fx((t+(delt/2)),(x+(delt*vx)/2),(y+(delt*vy)/2),(z+(delt*vz)/2),(vx+(k1x/2)),(vy+(k1y/2)),(vz+(k1z/2)));
k2y = delt*fy((t+(delt/2)),(x+(delt*vx)/2),(y+(delt*vy)/2),(z+(delt*vz)/2),(vx+(k1x/2)),(vy+(k1y/2)),(vz+(k1z/2)));
k2z = delt*fz((t+(delt/2)),(x+(delt*vx)/2),(y+(delt*vy)/2),(z+(delt*vz)/2),(vx+(k1x/2)),(vy+(k1y/2)),(vz+(k1z/2)));

k2x

k3x = delt*fx((t+(delt/2)),(x+((delt*vx)/2)+((delt*k1x)/4)),(y+((delt*vy)/2)+((delt*k1y)/4)),(z+((delt*vz)/2)+((delt*k1z)/4)),(vx+(k2x/2)),(vy+(k2y/2)),(vz+(k2z/2)));
k3y = delt*fy((t+(delt/2)),(x+((delt*vx)/2)+((delt*k1x)/4)),(y+((delt*vy)/2)+((delt*k1y)/4)),(z+((delt*vz)/2)+((delt*k1z)/4)),(vx+(k2x/2)),(vy+(k2y/2)),(vz+(k2z/2)));
k3z = delt*fz((t+(delt/2)),(x+((delt*vx)/2)+((delt*k1x)/4)),(y+((delt*vy)/2)+((delt*k1y)/4)),(z+((delt*vz)/2)+((delt*k1z)/4)),(vx+(k2x/2)),(vy+(k2y/2)),(vz+(k2z/2)));

k4x = delt*fx((t+delt),(x+(delt*vx)+((delt*k2x)/2)),(y+(delt*vy)+((delt*k2y)/2)),(z+(delt*vz)+((delt*k2z)/2)),(vx+k3x),(vy+k3y),(vz+k3z));
k4y = delt*fy((t+delt),(x+(delt*vx)+((delt*k2x)/2)),(y+(delt*vy)+((delt*k2y)/2)),(z+(delt*vz)+((delt*k2z)/2)),(vx+k3x),(vy+k3y),(vz+k3z));
k4z = delt*fz((t+delt),(x+(delt*vx)+((delt*k2x)/2)),(y+(delt*vy)+((delt*k2y)/2)),(z+(delt*vz)+((delt*k2z)/2)),(vx+k3x),(vy+k3y),(vz+k3z));

%updating time
t = t + delt;

%updating position
x(index+1) = x(index) + (delt*vx(index)) + ((delt*(k1x + k2x + k3x))/6);
y(index+1) = y(index) + (delt*vy(index)) + ((delt*(k1y + k2y + k3y))/6);
z(index+1) = z(index) + (delt*vz(index)) + ((delt*(k1z + k2z + k3z))/6);

%updating velocity
vx(index+1) = vx(i) + ((k1x + 2*k2x + 2*k3x + k4x)/6);
vy(index+1) = vy(i) + ((k1y + 2*k2y + 2*k3y + k4y)/6);
vz(index+1) = vz(i) + ((k1z + 2*k2z + 2*k3z + k4z)/6);

index = index + 1;
end

vx
vy
vz

x
y
z