How to create a flow chart?

How would I make a flow chart that describes this code?

clear all; close all; clc;

%inputs
d1f = 60; %diameter of the blades in feet 
r1 = (d1f/2)/3.281; %radius of the blades in meters
Vw1 = 24.51; %average wind speed in miles per hour
Vw = Vw1/2.237; %converting the wind speed to meters per second
%defining other constants
rho = 1.2754; %density of air in kg/m^3
Cq = 0.27; %estimates using diameter and torque from systems project
T1 = (1/2)*rho*pi*(r1^3)*Cq*(Vw^2);
V1 = 89;
%defining needed values
tspan = [0, 25]; % defining 
x0 = [0;0;0;0;0]; %defining initial conditions

[t, x] = ode15s(@(t,x) odefun(t,x,T1,V1),tspan,x0) %using ode15s to solve


u = [T1;V1]; %defines the u matrix
%defining variables
m = 18130.59;  %equivalent mass including blades, rotor and low speed shaft in (kg)
k_t = 80363.83655; %stiffness of the low speed shaft, found from textbook¹ (N/m)
c_aero = 4561000.755; %found from doing (T1/V1)*400 (Ns)
c_t = 4561000.755; %damping of the motor, estimated from doing multiple runs (Ns)
k_b = 6.86;  %stiffness of the motor³ (N/m)
R = 0.48;%10.41; %resistance found from motor spec sheet⁴ (Ohms)
L = 3.4;%0.644;  %thermal inductance found from spec sheet⁴ (Hertz)
I1 =  134000;  %moment of inertia of rotor/blade assembly, estimated (kg*m^2)
I2 =  0.000027525; %moment of inertia of the motor, found using mr^2 (kg/m^2)
N = 5.71;   %gear ratio⁵
e= 2.1e11;  %young's modulus of the material (Pa)
d_l = 3;%0.7736; %diameter of low speed shaft (meters)⁶
d_h = 0.3;%0.007;  %diameter of the high speed shaft (meters)⁸
l_l = 6;%3.5; %length of the low speed shaft (meters)⁶
l_h=0.6;%0.3;  %length of the high speed shaft (meters)⁸
area_l = pi*(d_l/2)^2; %area of the low speed shaft (m^2)
area_h = pi*(d_h/2)^2;  %area of the high speed shaft (m^2)
k_l = ((e*area_l)/l_l); %stiffness of low speed shaft (N/m)
k_h = ((e*area_h)/l_h); %stiffness of the high speed shaft (N/m)
%State matrices
A = [0 1 0 0 0
  (-N*k_h)/(I1) -c_aero/(I1) k_h/(I1) 0 0
  0 0 0 1 0
  (N*k_h)/(I2) 0 -k_h/I2 -c_t/I2 k_t/I2
  0 0 0 -k_b/L -R/L];
B = [0 0
  1/(I1) 0
   0 0
   0 0
   0 -1/L];
C = eye(5);
D = [0];

% For the fifth row (current), compute the 
for i = 1:length(x)
dxdt(:,i) = A*x(i,:)'+B*u;%sets up the state matrix equation
end
el = dxdt(5,:);
P = el*V1;
w1 = integral(@(el) el.*V1,0,inf)
clear m k_t c_aero c_t k_b R L I1 I2 N e d_l d_h l_l l_h area_l area_h k_l k_h C D

figure()
plot(t,P)
title('Power vs. Time')
figure()
plot(t,w1,'-o')
title('Work vs. Time')

figure
%plots
plot(t,x)
grid
xlabel('Time')
ylabel('Amplitude')
NrSp = size(x,2);
figure
%for loop creating an array of sub plots for x1 to x5
for k = 1:NrSp  %NrSp is the number of subplots equal to the column size of x
    %plotting code
  subplot(NrSp,1,k)
  plot(t, x(:,k))
  grid
  title(sprintf('x_%d',k))
end
xlabel('Time')

myfun = @(t,x) A*[x(1);x(2);x(3);x(4);x(5)]+B*u
 




function [dxdt] = odefun (t,x,T1,V1) %sets up a function using the state matrices
u = [T1;V1]; %defines the u matrix
%defining variables
m = 18130.59;  %equivalent mass including blades, rotor and low speed shaft in (kg)
k_t = 80363.83655; %stiffness of the low speed shaft, found from textbook¹ (N/m)
c_aero = 4561000.755; %found from doing (T1/V1)*400 (Ns)
c_t = 4561000.755; %damping of the motor, estimated from doing multiple runs (Ns)
k_b = 6.86;  %stiffness of the motor³ (N/m)
R = 0.48;%10.41; %resistance found from motor spec sheet⁴ (Ohms)
L = 3.4;%0.644;  %thermal inductance found from spec sheet⁴ (Hertz)
I1 =  134000;  %moment of inertia of rotor/blade assembly, estimated (kg*m^2)
I2 =  0.000027525; %moment of inertia of the motor, found using mr^2 (kg/m^2)
N = 5.71;   %gear ratio⁵
e= 2.1e11;  %young's modulus of the material (Pa)
d_l = 3;%0.7736; %diameter of low speed shaft (meters)⁶
d_h = 0.3;%0.007;  %diameter of the high speed shaft (meters)⁸
l_l = 6;%3.5; %length of the low speed shaft (meters)⁶
l_h=0.6;%0.3;  %length of the high speed shaft (meters)⁸
area_l = pi*(d_l/2)^2; %area of the low speed shaft (m^2)
area_h = pi*(d_h/2)^2;  %area of the high speed shaft (m^2)
k_l = ((e*area_l)/l_l); %stiffness of low speed shaft (N/m)
k_h = ((e*area_h)/l_h); %stiffness of the high speed shaft (N/m)
%State matrices
A = [0 1 0 0 0
  (-N*k_h)/(I1) -c_aero/(I1) k_h/(I1) 0 0
  0 0 0 1 0
  (N*k_h)/(I2) 0 -k_h/I2 -c_t/I2 k_t/I2
  0 0 0 -k_b/L -R/L];
B = [0 0
  1/(I1) 0
   0 0
   0 0
   0 -1/L];
C = eye(5);
D = [0];

dxdt = A*x+B*u;%sets up the state matrix equation

end