I writing program for rung kutta th4 now i have a problem who can help me?

% clear;
clc;
%%         precision of work
    n=100;                                   %steps 
    dr=1/n;                                   %minimum step radius
    r=1/n:dr:1+1/n;                       %0=<r=<1 radius
    %%     
    vdc=2.4;                                  %max voltage(v)
    dv=0.001;                                %step voltage(v)
    %%     capasitor properties
    R=500e-6;                                %max radius plate(m)
    h=1e-6;                                     %thickness plate(m)
    g=2e-6;                                     %Initial gap capasitor(m)
    %%     material properties
    E=150e9;                                  %Young’s modulus (pa)
    e0=8.854e-12;                          %Electrical permittivity of air
    nu=0.06;                                    %Poisson ratio
    rho=2300;                                  %Density (kg/m^3)    
    %%     Fixed functions
    D=(E*h^3)/(12*(1-nu^2));           %flexural strength
    %tstar=R^2*sqrt((rho*h)/D);
    %betha2=(c*R^4)/(D/tstar);
    alpha=(e0*R^4)/(2*D*g^3); 
    %%     mode shape
    phi=cos((pi*r)/2).^2;
    phi1=-cos(pi.*r/2).*pi.*sin(pi.*r/2);
    phi2=((pi^2*sin(pi*r/2).^2)/2)-((cos(pi.*r/2).^2*pi^2)/2);
    phi3=pi^3*sin(pi.*r/2).*cos(pi.*r/2);
    phi4=((pi^4*cos(pi*r/2).^2)/2)-((pi^4*sin(pi*r/2).^2)/2);
    %% 
    ws=zeros(1,n+1);                %zeros matrix for ws
    wprim=zeros(1,n+1);
    V=0;                                     %first voltage
    zita=0.01;                            %Damping ratio=C/C(critical)  %C(critical)=2*sqrt(k*m)
    %% time option
    T=1;
    dt=0.01;       %step time
    t=0;          %first time
    y1=zeros(1,(T/dt)+1); y1=0;
    y2=zeros(1,(T/dt)+1); y2=0;   
    %%     loop
    for i=1:vdc/dv
    
        %%  static
    Tstr=((E*h)/2*R*(1-nu^2)).*sum((wprim.^2)*dr);        %stretching force
    betha1=Tstr.*(g^2*R)/D;
    
    Fs=sum(((2*alpha.*V.*dv.*phi)./(1-ws).^2)*dr);
    kelec=sum(((2*alpha.*V.^2.*phi.^2)./((1-ws).^3))*dr);
    kmech=sum(((phi).*(phi4+(2./r).*phi3-((1./r.^2).*phi2)+((1./r.^3).*phi1)+...
        (betha1.*(phi2+((1./r).*phi1))))).*dr);
    keq=kmech-kelec;
    as=Fs/keq;
    psi=as*phi;
    ws=ws+psi;
    
    for j=1:n
        wprim=(ws(j+1)-ws(j))/dr;     %dw
    end
    
    V=V+dv;
wss(i)=ws(1);
VV(i)=V;

    %%  dynamic
    m=sum((phi.^2)*dr);
    k=sum(((phi).*(phi4+(2./r).*phi3-((1./r.^2).*phi2)+((1./r.^3).*phi1)+...
        (betha1.*(phi2+((1./r).*phi1))))).*dr);
    c=2*zita*sqrt(k*m);
    end
    %% runge kutta loop
      for a=1:T/dt
          %  first
          F=sum(((alpha.*(V.^2).*phi)./(1-(y1.*phi)).^2)*dr);
          s1=y2;
          p1=(1/m)*(F-(c*y2)-(k*y1));
             %  second
             F=sum(((alpha.*(V.^2).*phi)./(1-((y1+(dt/2)*s1).*phi)).^2)*dr);
             s2=y2+(dt/2)*p1;
             p2=(1/m)*(F-(c*(y2+(dt/2)*p1))-(k*(y1+(dt/2)*s1)));
                %  third
                F=sum(((alpha.*(V.^2).*phi)./(1-((y1+(dt/2)*s2).*phi)).^2)*dr);
                s3=y2+(dt/2)*p2;
                p3=(1/m)*(F-(c*(y2+(dt/2)*p2))-(k*(y1+(dt/2)*s2)));
                   % fourth
                   F=sum(((alpha.*(V.^2).*phi)./(1-((y1+(dt*s3)).*phi)).^2)*dr);
                   s4=y2+(dt*p3);
                   p4=(1/m)*(F-(c*(y2+(dt*p3)))-(k*(y1+(dt*s3))));
                   
          y1=y1+(dt/6)*(s1+(2*s2)+(2*s3)+s4);
          y2=y2+(dt/6)*(p1+(2*p2)+(2*p3)+p4);
          
          t=t+dt;
          tt(a)=t;
          y1(a)=y1;
          y2(a)=y2;
      end
    
    %% 
    

    
%plot(tt,y1.*phi)
%plot(tt,y2.*phi)
plot(y1.*phi,y2.*phi)
%plot(VV,1-wss)
grid on

[Walter Roberson]

    for j=1:n
        wprim=(ws(j+1)-ws(j))/dr;     %dw
    end

That overwrites all of wprim each iteration of the loop.