// eqondes.sce
//
clear all
n = 200;
A = 2*diag(ones(n,1))-diag(ones(n-1,1), -1) ...
    -diag(ones(n-1,1), 1);
A = -(n+1)^2 * A;
[R, D] = spec(A);
D = diag(D);

function y=g(x)
    y = 27/4 * x.^2 .* (1-x);
endfunction

h = 1.0 / (n+1);
x = (h : h : 1-h)'; // xj
u0 = g(x);
plot(x, u0, '.-'); xgrid();
dt = 0.02;
t = 0; T =15;
//
while (t<T)
    t = t + dt;
    u = 0*x;
    for j=1:n
        u = u + (u0'*R(:,j)) * cos(sqrt(-D(j))*t) * R(:,j);
    end; // for j    
    drawlater;
    clf; plot(x, u, '.-', 0, -1.2, '.', 1, 1.2, '.'); xgrid();
    drawnow
    disp(t)
end // while