function m = supersize
    global e k p

    e = 5000;
    k = 7800;
    p = 21.6;

    m0 = 84;
    tspan = [0,100];
    opts = odeset('RelTol',1.0e-1);
    [t,m] = ode45(@mbilanz,tspan,m0,opts);

    mexact = (e/p) + (m0 - e/p)*exp(-t*p/k);

    subplot(1,2,1)
    plot(t,mexact,'b')
    subplot(1,2,2)
    plot(t,m,'r')

end

function dm = mbilanz(t,m)

    global e k p

    dm = (e/k) - (p/k)*m;

end
