How can I speed up my code? I've looked through the suggestions, but I feel like none of them apply.

I'm just doing a bunch of calculations. What I have is a Classic Runge-Kutta method of solving a system of second order differential equations. Each iteration of the loop has to do many calculations. I need to run this code at a very fine dt with hundreds of thousands of iterations, which would take approximately 2 days by my estimation. I'm pretty new to matlab, so I feel like my code is super clunky and inefficient, what can I do to improve this?

function classicRK(dt,n,Q,a0)
syms h2 a2 h1 a1 h a %Define variables. In this notation, the numbers represent the order of the derivative of that variable with respect to time.
M_hh=1; %Define the parameters for the wing.
M_ha=0.625;
M_aa=1.25;
M_ah=0.25;
D_h=0.1;
D_a=0.25;
K_h=0.2;
K_a=1.25;
k_NL=10;
L=Q*a;
M=-0.7*Q*a;


%These are the equations that govern the motion of the wing.
eqn1=(M_hh*h2)+(M_ha*a2)+(D_h*h1)+(K_h*h)+L==0;
eqn2=(M_aa*a2)+(M_ah*h2)+(D_a*a1)+(K_a*(1+(k_NL*h^2))*a)+M==0;

%Here, the capital letters H and A are substitutions to get rid of the 2nd order derivatives.

[H1(a1,a,h,h1),A1(a,a1,h,h1)]=solve(eqn1,eqn2,h2,a2); %solves the system for the 2nd order derivatives.
%H1(1,2,3,4)
%A1(1,2,3,4)

t_val=zeros(1,n); %initial values
h_val=zeros(1,n); 
H_val=zeros(1,n);
A_val=zeros(1,n);
a_val=zeros(1,n);

a_val(1)=a0; %0<a<0.08 radians


for i=2:n
    t_val(i)=t_val(i-1)+dt;
    
    if(mod(i,100)==0) %tracker
        disp(i);
    end
    
    ka1=dt*A_val(i-1);
    kh1=dt*H_val(i-1);
    ka2=(1/2)*ka1;
    kh2=(1/2)*kh1;
    ka3=(1/2)*ka2;
    kh3=(1/2)*kh2;
    ka4=ka3+ka1;
    kh4=kh3+kh1;
    
    kA1=dt*A1(a_val(i-1),A_val(i-1),h_val(i-1),H_val(i-1));
    kH1=dt*H1(A_val(i-1),a_val(i-1),h_val(i-1),H_val(i-1));
    kA2=dt*A1(a_val(i-1)+ka2,A_val(i-1)+(kA1*(1/2)),h_val(i-1)+kh2,H_val(i-1)+(kH1*(1/2)));
    kH2=dt*H1(A_val(i-1)+(kA1*(1/2)),a_val(i-1)+ka2,h_val(i-1)+kh2,H_val(i-1)+(kH1*(1/2)));
    kA3=dt*A1(a_val(i-1)+ka3,A_val(i-1)+(kA2*(1/2)),h_val(i-1)+kh3,H_val(i-1)+(kH2*(1/2)));
    kH3=dt*H1(A_val(i-1)+(kA2*(1/2)),a_val(i-1)+ka3,h_val(i-1)+kh3,H_val(i-1)+(kH2*(1/2)));
    kA4=dt*A1(a_val(i-1)+ka4,A_val(i-1)+(kA3),h_val(i-1)+kh4,H_val(i-1)+(kH3));
    kH4=dt*H1(A_val(i-1)+(kA3),a_val(i-1)+ka4,h_val(i-1)+kh4,H_val(i-1)+(kH3));
    
    H_val(i)=H_val(i-1)+((1/6)*kH1)+((1/3)*kH2)+((1/3)*kH3)+((1/6)*kH4);
    A_val(i)=A_val(i-1)+((1/6)*kA1)+((1/3)*kA2)+((1/3)*kA3)+((1/6)*kA4);
    a_val(i)=a_val(i-1)+((1/6)*ka1)+((1/3)*ka2)+((1/3)*ka3)+((1/6)*ka4);
    h_val(i)=h_val(i-1)+((1/6)*kh1)+((1/3)*kh2)+((1/3)*kh3)+((1/6)*kh4);    

end

writeData('classicRK.csv',t_val,a_val,h_val)
plot(t_val,a_val)
ylim([-0.1,0.1])
xlim([0,60])

[dpb]

Why "roll your own" when MATLAB supplies builtin toolsets? <Ordinary-differential-equations>

[Walter Roberson]

I would suggest to you that using variables a_val and also A_val is more difficult to read than it needs to be; and likewise for other variables that you have used that uppercase / lowercase pattern with. Remember that other people need to read your code.

There are a number of re-used expressions in your code, such as a_val(i-1): you should assign those to variables so that MATLAB does not need to recalculate them each place they occur.