function ks,


    sstr_v=[];
    for j=1:16,
        if j<10,
            sstr_v(j,:)=['vQ','0',num2str(j)];  
        else
            sstr_v(j,:)=['vQ',num2str(j)];
        end
    end
    
    sstr_Q=[];
    for j=1:16,
        if j<10,
            sstr_Q(j,:)=['Q','0',num2str(j)];
        else
            sstr_Q(j,:)=['Q',num2str(j)];
        end
    end
    

    sstr_kQ1=[];
    for j=1:16,
        sstr_kQ1(j,:)=['k',sstr_Q(j,:),'1'];
    end
        sstr_kQ2=[];
    for j=1:16,
        sstr_kQ2(j,:)=['k',sstr_Q(j,:),'2'];
    end
        sstr_kQ3=[];
    for j=1:16,
        sstr_kQ3(j,:)=['k',sstr_Q(j,:),'3'];
    end
        sstr_kQ4=[];
    for j=1:16,
        sstr_kQ4(j,:)=['k',sstr_Q(j,:),'4'];
    end
   
    for j=1:16,
        sstr_Revised_Q(j,:)=[sstr_Q(j,:),'(i)'];
    end    
    str_Revised_Q='Q00(i)';
    for j=1:16,
        str_Revised_Q=[str_Revised_Q,',',sstr_Revised_Q(j,:)];
    end    
    
    for j=1:16,
        sstr_Revised_Q_1(j,:)=[sstr_Q(j,:),'(i)+',sstr_kQ1(j,:),'/2'];
    end    
    str_Revised_Q_1='Q00(i)+kQ001/2';
    for j=1:16,
        str_Revised_Q_1=[str_Revised_Q_1,',',sstr_Revised_Q_1(j,:)];
    end
    
    for j=1:16,
        sstr_Revised_Q_2(j,:)=[sstr_Q(j,:),'(i)+',sstr_kQ2(j,:),'/2'];
    end        
    str_Revised_Q_2='Q00(i)+kQ002/2';
    for j=1:16,
        str_Revised_Q_2=[str_Revised_Q_2,',',sstr_Revised_Q_2(j,:)];
    end    
    
    for j=1:16,
        sstr_Revised_Q_3(j,:)=[sstr_Q(j,:),'(i)+',sstr_kQ3(j,:)];
    end        
    str_Revised_Q_3='Q00(i)+kQ003';
    for j=1:16,
        str_Revised_Q_3=[str_Revised_Q_3,',',sstr_Revised_Q_3(j,:)];
    end     
    
    str_v='vQ00';
    for j=1:16,
        str_v=[str_v,',',sstr_v(j,:)];
    end
    str_Q='Q00';
    for j=1:16,
        str_Q=[str_Q,',',sstr_Q(j,:)];
    end

    Q00=1;
    for j=1:16,
    eval([sstr_Q(j,:),'(1)','=1;']);   %initializing
    end
    plot(Q01,'b*');
    Q01
    hold on;
    dt=0.0001;
    i=0;
    mu=4;
while(i<10),
    i=i+1;
    eval(['[',str_v,']=vel(mu,',str_Revised_Q,');']); %calculate v's the 1st time
    kQ001=dt*vQ00;
    for j=1:16,
        eval([sstr_kQ1(j,:),'=dt*',sstr_v(j,:),';'])  %update kQ??1
    end
    eval(['[',str_v,']=vel(mu,',str_Revised_Q_1,');']); %calculate v's the 2nd time
    kQ002=dt*vQ00;    
    for j=1:16,
        eval([sstr_kQ2(j,:),'=dt*',sstr_v(j,:),';'])  %update kQ??2
    end
    eval(['[',str_v,']=vel(mu,',str_Revised_Q_2,');']); %calculate v's the 3rd time
    kQ003=dt*vQ00;        
    for j=1:16,
        eval([sstr_kQ3(j,:),'=dt*',sstr_v(j,:),';'])  %update kQ??3
    end
    eval(['[',str_v,']=vel(mu,',str_Revised_Q_3,');']); %calculate v's the 4th time
    kQ004=dt*vQ00;        
    for j=1:16,
        eval([sstr_kQ4(j,:),'=dt*',sstr_v(j,:),';'])  %update kQ??4
    end    
    %update Q's for the i'th time
    Q00(i+1)=Q00(i)+kQ001/6+kQ002/3+kQ003/3+kQ004/6;
    for j=1:16,
        eval([sstr_Q(j,:),'(i+1)=',sstr_Q(j,:),'(i)+',sstr_kQ1(j,:),'/6+',sstr_kQ2(j,:),'/3+',sstr_kQ3(j,:),'/3+',sstr_kQ4(j,:),'/6;']);
    end
    plot(Q01,'b*');
end




function [vQ0,vQ1,vQ2,vQ3,vQ4,vQ5,vQ6,vQ7,vQ8,vQ9,vQ10,vQ11,vQ12,vQ13,vQ14,vQ15,vQ16]=vel(mu,Q0,Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,Q10,Q11,Q12,Q13,Q14,Q15,Q16),
for k=0:16,
    sig=0;
    for i=0:k ,
        eval(['sig=sig+Q',num2str(i),'*','Q',num2str(k-i),';']);
    end
    if k~=16,
        for i=(k+1):16,
            eval(['sig=sig-2*Q',num2str(i),'*','Q',num2str(16-i),';']);
        end
    end
    eval(['vQ',num2str(k),'=(k*k-mu*k*k*k*k)*Q',num2str(k),'-k*sig;']);
end