####THIS PROGRAM RUNS IN SAGEMATH 8.1

def graph_mup1(p,q,r,N):
    import itertools
    K=["".join(seq) for seq in itertools.product("012",repeat=N)]
    T=[(x/5).function(x),(3*x/5+1/5).function(x),(x/5+4/5).function(x)]
    R=[0]
    W=[]
    mup=0
    for i in range(3^N):
        aux0=0
        aux=1
        paux=1
        for j in range(N):
            aux0=T[int(K[i][N-j-1])](aux0)
            aux=T[int(K[i][N-j-1])](aux)
            if K[i][N-j-1]=='0':
                paux=paux*p
            if K[i][N-j-1]=='1':
                paux=paux*q
            if K[i][N-j-1]=='2':
                paux=paux*r
        mup=mup+paux
        R.append(mup)        
        if i<3^N-1:
            R.append(mup)
        W.append(aux0)             
        W.append(aux)
    x_coords = W
    y_coords = R
    result=list_plot(zip(x_coords, y_coords),plotjoined=True,color='green')
    return result

def graph_mup2(p,q,r,N):
    import itertools
    K=["".join(seq) for seq in itertools.product("012",repeat=N)]
    T=[(x/5).function(x),(3*x/5+1/5).function(x),(x/5+4/5).function(x)]
    R=[0]
    W=[]
    mup=0
    for i in range(3^N):
        aux0=0
        aux=1
        paux=1
        for j in range(N):
            aux0=T[int(K[i][N-j-1])](aux0)
            aux=T[int(K[i][N-j-1])](aux)
            if K[i][N-j-1]=='0':
                paux=paux*p
            if K[i][N-j-1]=='1':
                paux=paux*q
            if K[i][N-j-1]=='2':
                paux=paux*r
        mup=mup+paux
        R.append(mup)        
        if i<3^N-1:
            R.append(mup)
        W.append(aux0)             
        W.append(aux)
    x_coords = W
    y_coords = R
    result=list_plot(zip(x_coords, y_coords),plotjoined=True,color='red')
    return result




f1(x)=1/5*x
f2(x)=3/5*x+1/5
f3(x)=1/5*x+4/5
plot(f1,0,1)+plot(f2,0,1)+plot(f3,0,1)
plot(f1.diff(),0,1,color='red')+plot(f2.diff(),0,1,color='green')+plot(f3.diff(),0,1)

a=graph_mup1(1/4,1/3,5/12,8)
b=graph_mup2(1/6,1/4,7/12,8)
p=a+b
p