clear all
nx = 3
ny = 6
%prompt = 'input nx'
%nx = input(prompt)
%prompt = 'input ny'
%ny = input(prompt)
nxint = 5
nyint = 5
ydomain = [0.0469100771 0.2307653450 0.5 0.7692346551 0.9530899230]
xdomain = 0.5*ydomain;
W = [0.1184634425 0.2393143353 0.2844444444 0.2393143353 0.1184634425]
p = nx*ny

% find the solution to the variational problem
for i = 1:nx
    for k = 1:nx
        for j = 1:ny
            for l = 1:ny
                p = (i-1)*ny + j;
                q = (k-1)*ny + l;
                sum1 = 0; 
                %ux^2
                temp = 1/(2*j+2*l-3)-2/(2*j+2*l-1)+1/(2*j+2*l+1);
                sum1 = sum1 - temp*0.25*((2*i-2)*2*k+(2*k-2)*2*i)*0.5^(2*i+2*k-3)/(2*i+2*k-3);
                sum1 = sum1 + temp*(2*i)*(2*k)*0.5^(2*i+2*k-1)/(2*i+2*k-1);
                if (i*k>1)
                    %if i*k=1 the following term is zero (not there) if i = 1 or k = 1
                    sum1 = sum1 + temp*0.0625*(2*i-2)*(2*k-2)*0.5^(2*i+2*k-5)/(2*i+2*k-5);
                end
                sum1;
                %uy^2
                sum2 = 0;
                temp = 0.0625*0.5^(2*i+2*k-3)/(2*i+2*k-3)-0.5*0.5^(2*i+2*k-1)/(2*i+2*k-1)+0.5^(2*i+2*k+1)/(2*i+2*k+1);
                sum2 = temp*2*j*2*l/(2*j+2*l-1);
                sum2 = sum2 - temp*((2*j-2)*2*l+(2*l-2)*2*j)/(2*j+2*l-3);
                if (j*l>1)
                    %if j=l=1 the following term is zero (not there) if j = 1 or l = 1
                    sum2 = sum2 + temp*(2*j-2)*(2*l-2)/(2*j+2*l-5);
                end
                sum1
                sum2
                A(p,q) = sum1 + sum2 
            end
        end
    end
end
%2u
for i = 1:nx
    for j = 1:ny
        p = (i-1)*ny + j;
        %u
        term1 = 0.25*0.5^(2*i-1)/(2*i-1)-0.5^(2*i+1)/(2*i+1);
        term2 = 1/(2*j-1)-1/(2*j+1);
        f(p) = term1*term2;
    end
end
soln = linsolve(A,f')
A;
f

%Calculate the variation integral for the solution just derived
I = 0;
for i = 1:nx
    for k = 1:nx
        for j = 1:ny
            for l = 1:ny
                p = (i-1)*ny + j;
                q = (k-1)*ny + l;
                I = I - soln(p)*soln(q)*A(p,q);
            end
        end
    end
end
temp = 0.;
for i=1:nx
    for j = 1:ny
        p = (i-1)*ny + j;
        temp = temp + 2*soln(p)*f(p);
    end
end
temp %integral of 2u
I=I+temp %variational integral
error = I-0.02858521

% Calculate the variational integral for the solution just derived using
% Gaussian integration
for ii = 1:nxint
    x = xdomain(ii);
    for jj = 1:nyint
        y = ydomain(jj);
        sum1 = 0; sum2 = 0; sum3 = 0; 
        w = (1/4-x*x)*(1-y*y); 
        wx = -2*x*(1-y*y);
        wy = (1/4-x*x)*(-2*y);
        for i = 1:nx
            for j = 1:ny
                index = (i-1)*ny+j;
                sum1 = sum1 + w*soln(index)*x^(2*i-2)*y^(2*j-2); %u
                sum2 = sum2 + wx*soln(index)*x^(2*i-2)*y^(2*j-2); %ux
                sum2 = sum2 + w*soln(index)*(2*i-2)*x^(2*i-1)*y^(2*j-2);
                sum3 = sum3 + wy*soln(index)*x^(2*i-2)*y^(2*j-2); %uy
                sum3 = sum3 + w*soln(index)*(2*j-2)*x^(2*i-2)*y^(2*j-1); 
            end
        end
        u(ii,jj) = sum1;
        ux(ii,jj) = sum2;
        uy(ii,jj) = sum3;
    end
end

sum1 = 0; sum2 = 0; sum3 = 0; sum4 = 0.;
% the factor of 0.5 is because the x Gauss points are spread from 0 to 0.5
% rather than 0 to 1
for ii = 1:nxint
    for jj = 1:nyint
        sum1 = sum1 + 0.5*W(ii)*W(jj)*2*u(ii,jj); % 2u
        sum2 = sum2 - 0.5*W(ii)*W(jj)*ux(ii,jj)*ux(ii,jj); % -ux^2
        sum3 = sum3 - 0.5*W(ii)*W(jj)*uy(ii,jj)*uy(ii,jj); % -uy^2
        sum4 = sum4 + 0.5*W(ii)*W(jj); %1
    end
end
sum1
%sum2
%sum3
sum4
soln(1)/4 %center value
I = (sum1 + sum2 + sum3)
      
