
function [F,q1,q2,p1,p2,C_cai]=DM_funcs_USM(t,R,k1)
     C_cai = 200;
%     [k1]=bindingconst(C_cai);
    k2 = 1.2387;
    g1 = 0.0756;
    gp1 = 0.0709;
    fp1 = 0.2838;
%     k1 = 0.9;
%     k2 = 0.1399;
%     fp1 = 14.4496;
%     gp1 = 3.6124;
%     g1 = 0.1340;


    r(1)=R(1);
    r(2)=R(2);
    r(3)=R(3);
    r(4)=R(4);
    r(5)=R(5);
    r(6)=R(6);
    r(7)=R(7);
    
    gam = 100;
%Mean for cdf
    p1=r(2)/r(1);
    p2=r(5)/r(4);

%Standard deviation for cdf
    q1=(sqrt((r(3)/r(1))-((r(2)/r(1)).^(2))));    
    q2=(sqrt((r(6)/r(4))-((r(5)/r(4)).^(2))));

%r,phi and I values for 1st PDE M1_lambda
    r0=-p1/q1;
    r1=(1-p1)/q1;

    phi0=0.5*(1+erf((r0-p1)/(q1*sqrt(2))));
    phi1=0.5*(1+erf((r1-p1)/(q1*sqrt(2))));

    I0=-(exp(-((-p1./q1).^(2))/2))/(sqrt(2*pi));
    I1=-(exp(-((((1-p1)./q1)).^(2))/2))/(sqrt(2*pi));
    
    %r,phi and I values for 2nd PDE M2_lambda
    r20=-p2/q2;
    r21=(1-p2)/q2;
    
    phi20=0.5*(1+erf((r20-p2)/(q2*sqrt(2))));
    phi21=0.5*(1+erf((r21-p2)/(q2*sqrt(2))));

    I20=-(exp(-((-p2./q2).^(2))/2))/(sqrt(2*pi));
    I21=-(exp(-((((1-p2)./q2)).^(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
    K0=phi20;
     K01=phi21;
%     K0inf=phi2inf;
    
    K10=((p2.*phi20)+(q2.*I20));
    K11=(p2.*phi21)+(q2.*I21);
%     K12=(p2.*phi2inf)+(q2.*I2in);
    
    K20=((p2.^(2)).*phi20)+((2*p2.*q2).*I20)+((q2.^(2)).*(phi20+(r20*I20)));
    K21=((p2.^(2)).*phi21)+((2*p2.*q2).*I21)+((q2.^(2)).*(phi21+(r21*I21)));
%     K22=((p2.^(2)))+((2*p2.*q2).*I2in)+((q2.^(2)));
    
    K30=(p2.^(3).*phi20)+((3.*p2.^(2).*q2).*I20)+((3.*p2.*q2.^(2)).*((phi20)+(r20*I20)))+((q2.^(3).*(2+r20.^(2)).*I20));
    K31=(p2.^(3).*phi21)+((3.*p2.^(2).*q2).*I21)+((3.*p2.*q2.^(2)).*(phi21+(r21*I21)))+(q2.^(3).*(2+(r21.^(2))).*I21);
%     K32=(p2.^(3))+((3.*p2.^(2).*q2).*I2in)+((3.*p2.*q2.^(2)));
    
    
% Components for the matrix F that will represent each moment,
%     M1_lambda and M2_lambda
  
    A0= ((fp1*(1-r(7)))/2)-(fp1*(J11-J10)*r(1));%-2*(fp1*(K11-K10)*r(4));
    A1=((fp1*(1-r(7)))/3)-(fp1*(J21-J20)*r(1));%-2*(fp1*(K21-K20)*r(4));
    A2=((fp1*(1-r(7)))/4)-(fp1*(J31-J30)*r(1));%-2*(fp1*(K31-K30)*r(4));
    
    B0=(3*(fp1+gp1)*J0)+gp1*(J11-J10)+(4*gp1)*(p1-J11);
    B1=(3*(fp1+gp1)*J10)+gp1*(J21-J20)+(4*gp1)*((p1.^(2)+q1.^(2))-J21);
    B2=(3*(fp1+gp1)*J20)+gp1*(J31-J30)+(4*gp1)*((p1.^(3)+3*p1*q1.^(2))-J31);
    
    C0=k1;
    C1=k1*p2;
    C2=k1*(p2.^(2)+q2.^(2));
    
    D0=k2;
    D1=k2*p1;
    D2=k2*(p1.^(2)+q1.^(2));
    
    E0=(20*g1*K0)+g1*(K11-K10)+(g1)*(1-K01);
    E1=(20*g1*K10)+g1*(K21-K20)+(g1)*((p2)-K11);
    E2=(20*g1*K20)+g1*(K31-K30)+(g1)*(p2.^(2)+q2.^(2)-K21);
    
    V = gam*(A1-E0-B0)/(1+gam*(D0+C0));
    
    F=[A0-B0*r(1)+C0*r(4)-k2*r(1);A1-B1*r(1)+C1*r(4)-k2*r(2)-V*r(1);A2-B2*r(1)+C2*r(4)-k2*r(3)-2*V*r(2);D0*r(1)-E0*r(4)-k1*r(4);D1*r(1)-E1*r(4)-k1*r(5)-V*r(4);D2*r(1)-E2*r(4)-k1*r(6)-2*V*r(5);-k1*r(7)+(1-r(7))*k2];
    