function [R,T,force,stiffness,em10,n1,P1,Q1,tt]=Zah(trange,N,M,xrange,f)
%[R,T,force,stiffness,em11,em10,em12,n1,P1,Q1,tt,vee]
%N number of time points
%trange specifies the interval of time steps
%T time grid
%X0 space grid
T=linspace(trange(1),trange(2),N);
X0=linspace(xrange(1),xrange(2),M);
%dt=timestep
dt=(trange(2)-trange(1))/(N);

%Initial conditions
M10=1;
M11=0.5;
M12=0.4;
R=[M10;M11;M12];
g=@(r) mis(r,N,M,f);%to pass into fsolve to find equilibria solutions
options = optimoptions('fsolve');
[R,~]=fsolve(g,R,options);

  
for j=1:N
    m10=R(1);
    m11=R(2);
    m12=R(3);
    
    %  t=0+j*dt;
    e=dt;
    tt(j)=0+j*e;
    
    
    X=X0;%-lengthvar(t,f);%update X for each L as L updates for each T
    
    % L=lengthvar(t,f);
    % L=(-25*sin(2*pi*f*t));
    % L=0;
    [alpha,V] =lengthvar(tt(j),f);
    vee(j)=V;
    fp1=0.06;
    g1=0.10;
    g2=0.0040;
    g3=0.01;
    
    %Mean for cdf
    p1=m11/m10;
    
    P1(j)=m11/m10;
    %p1=0.5;
    %p2=M21(i)/M20(i);
    %Standard deviation for cdf
    q1=sqrt(((m12/m10)-((m11/m10).^(2))));
    Q1(j)=sqrt(((m12/m10)-((m11/m10).^(2))));
    
    r0=-p1./q1;
    % r0=0;
    r1=(1-p1)./q1;
    %r1=1;
    rinf=(Inf-p1)./q1;
    %rinf=Inf;
    phi0=cdf('Normal',r0,p1,q1);
    % phi0=erf(-p1./q1);
    phi1=cdf('Normal',r1,p1,q1);
    %phi1=erf((1-p1)./q1);
    
    phinf=cdf('Normal',rinf,p1,q1);
    
    %phinf=erf((Inf-p1)./q1);
    
    I0=-(exp(-((-p1./q1).^(2))/2))/(sqrt(2*pi));
    I1=-(exp(-((((1-p1)./q1)).^(2))/2))/(sqrt(2*pi));
    I2=-(exp(-(((Inf-p1)./q1)).^(2))/2)/(sqrt(2*pi));
    
    
    
    %Functions for the rhs of the first PDE M1_lambda
    J0=phi0;
    J01=phi1;
    J0inf=phinf;
    
    J10=((p1.*phi0)+(q1.*I0));
    J11=(p1.*phi1)+(q1.*I1);
    J12=(p1.*phinf)+(q1.*I2);
    
    J20=((p1.^(2)).*phi0)+((2*p1.*q1).*I0)+((q1.^(2)).*(phi0+(r0*I0)));
    J21=((p1.^(2)).*phi1)+((2*p1.*q1).*I1)+((q1.^(2)).*(phi1+(r1*I1)));
    J22=((p1.^(2)))+((2*p1.*q1).*I2)+((q1.^(2)));
    
    J30=(p1.^(3).*phi0)+((3.*p1.^(2).*q1).*I0)+((3.*p1.*q1.^(2)).*((phi0)+(r0*I0)))+((q1.^(3).*(2+(r0.^(2))).*I0));
    J31=(p1.^(3).*phi1)+((3.*p1.^(2).*q1).*I1)+((3.*p1.*q1.^(2)).*((phi1)+(r1*I1)))+(q1.^(3).*(2+(r1.^(2))).*I1);
    J32=(p1.^(3))+((3.*p1.^(2).*q1).*I2)+((3.*p1.*q1.^(2)));
    
    %Functions defined for the RHS of the second PDE M2_lambda
    
    %Components for the matrix F that will represent each moment,
    %M1_lambda and M2_lambda
    
    %
    b0=fp1/2;
    b1=fp1/3;
    b2=fp1/4;
    F0=((g2*J0)+((fp1+g1)*(J11-J10))+(g1*(p1-J11))+(g3*(p1-J11-1+J01)))*m10;
    F1=((g2*J10)+((fp1+g1)*(J21-J20))+(g1*(p1.^(2)+q1.^(2)-J21))+(g3*((p1.^(2)+q1.^(2))-J21-p1+J11)))*m10;
    F2=((g2*J20)+((fp1+g1)*(J31-J30))+(g1*(p1.^(3)+3*p1*q1.^(2)-J31))+(g3*((p1.^(3)+3*p1*q1.^(2))-J31-(p1.^(2)+q1.^(2))+J21)))*m10;
    
    % Putting together the matrix F for the RHS of the system of
    % equations
    
    
    % F=[b0-F0*m10;b1-F1*m10;b2-F2*m10];
    F=[b0-F0;b1-F1-V*m10;b2-F2-2*V*m11];
    
    Rnew=R+dt*F;
    R=Rnew;
    %m10=r(1);
    
%  m=1;
    
    n1=(R(1)/((sqrt(2*pi))*q1)).*exp(-((X-p1).^(2))./(2*(q1.^(2))));
    force(j)=trapz(X,X.*(n1));
    stiffness(j)=trapz(X,n1);
     
    em11(j)= R(2);
    em10(j)=R(1);
    em12(j)=R(3);
    
%    % m=ones(1,N/5000);
%     t=zeros(1,N/5000);
% %     em11=zeros(1,N/5000);
%      em10=zeros(1,N/5000);
% %     em12=zeros(1,N/5000);
% 
%     if (mod(j,5000)==0)% && j>=16500 && j<=16700) 
%        
%         em10(m)=R(2);
%        
% %           n=m+1;
% %      % t(m)=tt;   
% %     em11(1,n)= R(2);
% %     em10(n)=R(1);
% %     em12(n)=R(3);
% %         n=n+1;
%     end 
%      mnew=m+1;
%         m=mnew;
% %         
% %         figure(1)
% %         plot(X0,n1)
% %         drawnow
% %         hold on
% %        q1
% %        
% %         pause
% %         %     [
% %                     figure(2)
% %                     plot(T,R(1))
% %         
% %                     drawnow
% %                     hold on
% %           pause(2)
    
   
end

