Plotting the mean and standard derivation of many stair-step graphs

I want to find the mean and standard derivation of several stairstep graphs with different length of steps.

Here is my function with some random variables.

clear all

%set parameter values
am=1;
bm=1;
ap=1;
bp=1;
ohm=100;
%set number of reactions
Tend=4000;
%set initial conditions for # molecules and for time t 
X=zeros(Tend,2);
t=zeros(Tend,1);
X(1,:)=[0 0]; t(1)=0;

j=1;
while t(j)<10
    %calculate propensities
    p1=am*ohm;
    p2=bm*X(j,1);
    p3=ap*X(j,1);
    p4=bp*X(j,2);
    psum=p1+p2+p3+p4;

    %select next reaction
    mu=rand(1);
    z1=0;z2=0;z3=0;z4=0;
    if 0 <= mu && mu < p1/psum % M generated
       z1=1;
        else if p1/psum <= mu && mu  < (p1+p2)/psum % M degraded
                z2=1;
            else if (p1+p2)/psum <= mu && mu < (p1+p2+p3)/psum % P generated
                    z3=1;
                else % P degraded
                     z4=1;
                end
            end
    end

    %update state
    X(j+1,1) = X(j,1) + z1 -z2;
    X(j+1,2) = X(j,2) + z3 -z4;
    %update time
    t(j+1) = t(j) + log(1/rand(1))/psum;
    %update counter
    j=j+1;
end

subplot(2,1,1)
stairs(t(1:j), X(1:j,1), 'k', 'linewidth', 2)
xlim([0 10])
ylabel('number of molecules')
xlabel('time')
title('M')

subplot(2,1,2)  
stairs(t(1:j), X(1:j,2), 'g', 'linewidth', 2)
xlim([0 10])
ylabel('number of molecules')
xlabel('time')
title('P')

This is the result of the first trial, and

this is the one of the second trial.

You can see that as the graph differs, j and the length of steps also differs everytime the function is plotted.

So I cannot use stdshade function since it assumes that 'the observation' occurs with uniform steps of the time.

I tried to find the value of M and P at t=0.01*n to make a new matrix, and then use stdshade.

while k<=1000
    while t(j)<10
        %calculate propensities
        p1=am*ohm;
        p2=bm*X(j,1);
        p3=ap*X(j,1);
        p4=bp*X(j,2);
        psum=p1+p2+p3+p4;

        %select next reaction
        mu=rand(1);
        z1=0;z2=0;z3=0;z4=0;
        if 0 <= mu && mu < p1/psum % M generated
           z1=1;
            else if p1/psum <= mu && mu  < (p1+p2)/psum % M degraded
                    z2=1;
                else if (p1+p2)/psum <= mu && mu < (p1+p2+p3)/psum % P generated
                        z3=1;
                    else % P degraded
                         z4=1;
                    end
                end
        end

        %update state
        X(j+1,1) = X(j,1) + z1 -z2;
        X(j+1,2) = X(j,2) + z3 -z4;
        %update time
        t(j+1) = t(j) + log(1/rand(1))/psum;
        %update counter
        j=j+1;
    end
    f1=stairs(t(1:j),X(1:j,1));
    f2=stairs(t(1:j),X(1:j,2));
    XData1=get(get(gca,'f1'),'XData');
    YData1=get(get(gca,'f1'),'YData');
    XData2=get(get(gca,'f2'),'XData');
    YData2=get(get(gca,'f2'),'YData');
    while n<=1000
        M1(k,n)=YData1(XData1==0.01*n);
        M2(k,n)=YData2(XData2==0.01*n);
        n=n+1;
    end
    n=1;
    j=1;
end
stdshade(M1,0:0.01:10)
stdshade(M2,0:0.01:10)

But it shows the error at 'matlab.graphics.axis.Axes/get'.

The expected result is as follows:

Can you help me? Thanks.

[dpb]

We can't run your function nor did you even show us what it produces -- not a lot to be done with such limited info.

That said, why can't you just use interp1 to create a set of uniform points for all the curves and then average those?

[kyuwon Oh]

Sorry for inconvinience. Now I attached the whole code and some figures. // I will try to use interp1. Thank you. // It worked! Thanks a lot for a brilliant advice.