%
% HW#1
%
% Amik St-Cyr NCAR 2007

clear all;

a=-pi;
b= pi;
N= 15; % Elements

x = a:(b-a)/(N-1):b;


[Operators] = Q1FEM(x);

%
% Solve the problem 
% 
% u - delta u = f
%
% Dirichlet u(-1)=u(1)=1
%

f = 1 + 2*sin(x);
f = Operators.M1DQ1*f';

A = Operators.M1DQ1 + Operators.K1DQ1;

% Impose boundary 
% conditions:
u    = zeros(N,1);
u(1) = 1;
u(N) = 1;
f    = f-A*u;

% Solve for interior homogeneous block:
u(2:N-1,1) = A(2:N-1,2:N-1)\f(2:N-1);

% Plot and compute error
plot(x,u)
norm(u - 1 - sin(x'))

%
% Do the same problem for Q2 Elements
% Using the same number of degrees 
% of freedom (points)
%
