Download<< Back
function run_bruss
N=100;
x = (1:N)/(N+1);
y0 = [1+sin((2*pi/(N+1))*(1:N)); repmat(3,1,N)];
t0 = 0;
T = 100;
options = odeset('Vectorized','on','JPattern',jpattern(N));
[t,y] = ode15s(@(t,x) bruss(t,x,N), [t0,T], y0, options);
uBC = ones(size(t));
vBC = 3*ones(size(t));
u = [uBC y(:,1:2:end) uBC];
v = [vBC y(:,2:2:end) vBC];
x = [0, x, 1];
mid = round(N/2);
figure;
plot(t, u(:,2), 'k', t, u(:,mid), 'b', t, u(:,N+1), 'r');
legend(['x = ' num2str(x(1))], ['x = ' num2str(x(mid))], ['x = ' num2str(x(N+1))]);
xlabel('t');
ylabel('u');
figure;
plot(t, v(:,2), 'k', t, v(:,mid), 'b', t, v(:,N+1), 'r');
legend(['x = ' num2str(x(1))], ['x = ' num2str(x(mid))], ['x = ' num2str(x(N+1))]);
xlabel('t');
ylabel('v');
figure;
surf(x, t, u);
xlabel('x');
ylabel('t');
zlabel('u');
figure;
surf(x, t, v);
xlabel('x');
ylabel('t');
zlabel('v');
function S = jpattern(N)
B = ones(2*N,5);
B(2:2:2*N,2) = zeros(N,1);
B(1:2:2*N-1,4) = zeros(N,1);
S = spdiags(B,-2:2,2*N,2*N);
end
end