function f = ctrl()

    clear

        syms kd kp k t x real syms ev es vd v0 sd real

            a = vd b = -1 / k*(v0 - vd)c = -k d = 1 / k*(v0 - vd) -
    sd

            solve(a * x + b * exp(c * x) + d)

        %
        solve(x * vd + (1 - exp(-k * x)) / k * (v0 - vd) - sd)

            syms x N real N = 2

    solve(0.5 * (1 + x ^ N) - 1 + (1 - x ^ N) / (1 - x), x)

        syms s(t) v(t)

            ode1 = diff(s) == v ode2 = diff(v) ==
    -kp*(v - vd)

        cond1 = s(0) == 0;
cond2 = v(0) == v0;

S = dsolve([ode1;ode2], [cond1;cond2])

sSol(t) = S.s
vSol(t) = S.v


syms v10 v20 v1d v2d k t a1 a2 real

syms a b c d ekt real

s1 = s_true(v10, v1d, k, t, ekt)
s2 = s_true(v20, v2d, k, t, ekt)
%s1 = v10*t + a1*t^2/2
%s2 = v20*t + a2*t^2/2

diff(s1,t)

e = simplify(a + s1^2 + s2^2 + b*s1 + c*s2 + d*s1*s2)
solve([e, ekt-exp(-k*t)], [t, ekt])


solve ( 
return

kp = 0.5;
k = 0.5;
vd = 10;
v0 = 5;

t = 0:0.1:4;

appr = 1-kp*t+(kp*t).^2/2;

%plot(t, t*vd + v0/kp - (exp(-kp*t).*(v0 - vd + vd*exp(kp*t)))/kp, t, exp(-kp*t).*(v0 - vd + vd*exp(kp*t)))
plot(t, t*vd + v0/kp - (exp(-kp*t).*(v0 - vd + vd*exp(kp*t)))/kp, t,t*vd + v0/kp - (appr.*(v0 - vd + vd./appr))/kp )

hold on



sd = 0:.1:30
t = (vd - v0 + k*sd + vd*lambertw(0, (exp(-(vd - v0 + k*sd)/vd)*(v0 - vd))/vd))/(k*vd);
plot(t, sd,'go')

%solve(t*vd + v0/kp - (exp(-kp*t)*(v0 - vd + vd*exp(kp*t)))/kp - sd, t)

ts = 0:0.1:4;
v0
vd
k

s = s_at_t(v0, vd, k, 4, 4)


ss = zeros(size(ts));
for i=2:length(ts),
    ss(i) = s_at_t(v0, vd, k, ts(i), 2)
end

plot(ts, ss)

legend('s','s_appr','s_lam','trap')

end


function s = s_true(v0, vd, k, t, ekt)
    
    s = t*vd + (v0 - vd)/k*(1 - ekt);
end

function v = v_at_t(v0, vd, kp, t)
    vd
    kp
    t
    v0
    v= vd + exp(-kp*t).*(v0 - vd);
end

function s = s_at_t(v0, vd, kp, t, n)
    dt = t/n;
    vs = v_at_t(v0, vd, kp, 0:dt:t)
    s = dt/2*sum(vs(1:n)+vs(2:n+1));
end