clear all;
clc;

z = linspace(0,2,100);   %Will be used for plotting functions

%Exact solution
syms y(x) x
Dy = diff(y);   %single derivative
D2y = diff(y,2); %double derivative
sol_ex(x) = dsolve(D2y + 3*Dy + cos(x)==0, y(0) == 0, Dy(2) == 1); % solving differential equation

for i = 1:100
    %y_ex(i) = -129.91*exp(-3*z(i))+129.81+0.1*cos(z(i))-0.3*sin(z(i)); % exact solution, in notes
    y_exact(i) = sol_ex(z(i));
end

plot(z,y_exact, 'ko');
ylabel('y(x)');
xlabel('x');
hold on;   %so that all graphs can be shown at the same time


% 2nd degree polynomial - only one unknown and only one trial function
syms x c1    % syms allows working with variables - in this case, they are x and c1.

N1(x) = x*(x-4);

y_star(x) = c1*N1(x)+x;   % writing approximiate solution

R(x) = diff(diff(y_star(x)))+ 3*diff(y_star(x)) + cos(x);  %Calculation of residual after putting y* in differential equation

% Putting it in Galerkin method with N1 = x(x-4)
F1 = int(N1(x)*R(x),0,2);   %intergrates between 0 and 2

%Solution for c1
[c1_1] = solve(F1);

%Plotting
for i = 1:100
    y_star2(i) = c1_1*N1(z(i))+z(i);
end
plot(z, y_star2, 'r');


%next guess - a 3rd order polynomial
syms x c1 c2
N1(x) = x*(x-4);
N2(x) = x*(x^2-12); 

y_star(x) = c1*N1(x)+c2*N2(x)+x;   % writing approximiate solution

R(x) = diff(diff(y_star(x)))+ 3*diff(y_star(x)) + cos(x);

% Putting it in Galerkin method with N1 = x(x-4) and N2 = x(x^2-12)
F1 = int(x*(x-4)*R(x),0,2);      % this is equal to (40*c1)/3 + (328*c2)/5 - 6*sin(2) - 12
F2 = int(x*(x^2-12)*R(x),0,2);   % this is equal to (232*c1)/5 + (1152*c2)/5 - 6*cos(2) - 28*sin(2) - 42 

%Solution for c1 and c2
[c1_2, c2_2] = solve(F1,F2);

%Plotting
for i = 1:100
    y_star3(i) = c1_2*N1(z(i))+c2_2*N2(z(i))+z(i);
end
plot(z, y_star3, 'b');


%next guess - a 4th order polynomial
syms x c1 c2 c3

N1(x) = x*(x-4);
N2(x) = x*(x^2-12); 
N3(x) = x*(x^3-32);

y_star(x) = c1*N1(x)+c2*N2(x)+c3*N3(x)+x;   % writing approximiate solution

R(x) = diff(diff(y_star(x)))+ 3*diff(y_star(x)) + cos(x);

% Putting it in Galerkin method with N1 = x(x-4), N2 = x(x^2-12) and N3 = x(x^3-32)
F1 = int(N1(x)*R(x),0,2);      
F2 = int(N2(x)*R(x),0,2);      
F3 = int(N3(x)*R(x),0,2);      

%Solution for c1, c2 and c3
[c1_3, c2_3, c3_3] = solve(F1,F2,F3); %somehow all values of ci's are negative here and hence this does not fit! Do I have bug in my computer?

%Plotting
for i = 1:100
    y_star4(i) = c1_3*N1(z(i))+c2_3*N2(z(i))+c3_3*N3(z(i))+z(i);    
end
%plot(z, y_star4, 'g');     %not plotted - something wrong 

%next guess - a 5th order polynomial
syms x c1 c2 c3 c4

N1(x) = x*(x-4);
N2(x) = x*(x^2-12); 
N3(x) = x*(x^3-32);
N4(x) = x*(x^4-80);

y_star(x) = c1*N1(x)+c2*N2(x)+c3*N3(x)+c4*N4(x)+x;   % writing approximiate solution

R(x) = diff(diff(y_star(x)))+ 3*diff(y_star(x)) + cos(x);

% Putting it in Galerkin method 
F1 = int(N1(x)*R(x),0,2);     
F2 = int(N2(x)*R(x),0,2);      
F3 = int(N3(x)*R(x),0,2);      
F4 = int(N4(x)*R(x),0,2);      

%Solution for ci's
[c1_4, c2_4, c3_4, c4_4] = solve(F1,F2,F3,F4);

%Plotting
for i = 1:100
    y_star5(i) = c1_4*N1(z(i))+c2_4*N2(z(i))+c3_4*N3(z(i))+c4_4*N4(z(i))+z(i);
end
plot(z, y_star5, 'k');

%next guess - a 6th order polynomial
syms x c1 c2 c3 c4 c5

N1(x) = x*(x-2*2);
N2(x) = x*(x^2-3*2^2); 
N3(x) = x*(x^3-4*2^3);
N4(x) = x*(x^4-5*2^4);
N5(x) = x*(x^5-6*2^5);

y_star(x) = c1*N1(x)+c2*N2(x)+c3*N3(x)+c4*N4(x)+c5*N5(x)+x;   % writing approximiate solution

R(x) = diff(diff(y_star(x)))+ 3*diff(y_star(x)) + cos(x);

% Putting it in Galerkin method
F1 = int(N1(x)*R(x),0,2);      
F2 = int(N2(x)*R(x),0,2);      
F3 = int(N3(x)*R(x),0,2);      
F4 = int(N4(x)*R(x),0,2);
F5 = int(N5(x)*R(x),0,2);

%Solution for ci's
[c1_5, c2_5, c3_5, c4_5, c5_5] = solve(F1,F2,F3,F4,F5);

%Plotting
for i = 1:100
    y_star6(i) = c1_5*N1(z(i))+c2_5*N2(z(i))+c3_5*N3(z(i))+c4_5*N4(z(i))+c5_5*N5(z(i))+z(i);  %something is messed up here too
end
%plot(z, y_star6, 'm');   %not plotted - something wrong 


%next guess - a 7th order polynomial
syms x c1 c2 c3 c4 c5 c6

N1(x) = x*(x-2*2);
N2(x) = x*(x^2-3*2^2); 
N3(x) = x*(x^3-4*2^3);
N4(x) = x*(x^4-5*2^4);
N5(x) = x*(x^5-6*2^5);
N6(x) = x*(x^6-7*2^6);

y_star(x) = c1*N1(x)+c2*N2(x)+c3*N3(x)+c4*N4(x)+c5*N5(x)+c6*N6(x)+x;   % writing approximiate solution

R(x) = diff(diff(y_star(x)))+ 3*diff(y_star(x)) + cos(x);

% Putting it in Galerkin method
F1 = int(N1(x)*R(x),0,2);      
F2 = int(N2(x)*R(x),0,2);      
F3 = int(N3(x)*R(x),0,2);      
F4 = int(N4(x)*R(x),0,2);
F5 = int(N5(x)*R(x),0,2);
F6 = int(N6(x)*R(x),0,2);

%Solution for ci's
[c1_6, c2_6, c3_6, c4_6, c5_6, c6_6] = solve(F1,F2,F3,F4,F5,F6);

%Plotting
for i = 1:100
    y_star7(i) = c1_6*N1(z(i))+c2_6*N2(z(i))+c3_6*N3(z(i))+c4_6*N4(z(i))+c5_6*N5(z(i))+c6_6*N6(z(i))+z(i);
end
plot(z, y_star7, 'c');

%next guess - a 8th order polynomial
syms x c1 c2 c3 c4 c5 c6 c7

N1(x) = x*(x-2*2);
N2(x) = x*(x^2-3*2^2); 
N3(x) = x*(x^3-4*2^3);
N4(x) = x*(x^4-5*2^4);
N5(x) = x*(x^5-6*2^5);
N6(x) = x*(x^6-7*2^6);
N7(x) = x*(x^7-8*2^7);

y_star(x) = c1*N1(x)+c2*N2(x)+c3*N3(x)+c4*N4(x)+c5*N5(x)+c6*N6(x)+c7*N7(x)+x;   % writing approximiate solution

R(x) = diff(diff(y_star(x)))+ 3*diff(y_star(x)) + cos(x);

% Putting it in Galerkin method
F1 = int(N1(x)*R(x),0,2);      
F2 = int(N2(x)*R(x),0,2);      
F3 = int(N3(x)*R(x),0,2);      
F4 = int(N4(x)*R(x),0,2);
F5 = int(N5(x)*R(x),0,2);
F6 = int(N6(x)*R(x),0,2);
F7 = int(N7(x)*R(x),0,2);

%Solution for ci's
[c1_7, c2_7, c3_7, c4_7, c5_7, c6_7, c7_7] = solve(F1,F2,F3,F4,F5,F6,F7);

%Plotting
for i = 1:100
    y_star8(i) = c1_7*N1(z(i))+c2_7*N2(z(i))+c3_7*N3(z(i))+c4_7*N4(z(i))+c5_7*N5(z(i))+c6_7*N6(z(i))+c7_7*N7(z(i))+z(i);
end
plot(z, y_star8, 'r--');

legend show;
legend('Location', 'Bestoutside');
legend BOXOFF;


