% ex3transmission1D To solve a transmission problem in 1D
%                   Pimentel & Sousa  05/25
%
close all
clc
clear all
% Define parameters
Nx = 200;                % Number of spatial points including boundaries
x0=-1;
xf=1;
L = xf-x0;               % Length of the domain
x = linspace(x0, xf, Nx);% Spatial grid
dx = L/(Nx - 1);         % Spatial step size
epsilon=dx;              % It should be larger than dx                       
theta1=2;
theta2=4;
% The F function
 F = @(p) p;
% Initial function for the Newton iteration
v0 = @(x) x/2;
plot(x,v0(x)','LineWidth',1.5);
hold on
u = v0(x)';
% Source term
for i=1:Nx;
    if x(i) > 0
       ff(i) = -x(i)^2; %troquei os sinais
    elseif x(i) < 0
       ff(i) = x(i)^4; %troquei os sinais
    end
end
% Exact solution
for i=1:Nx;
    if x(i) > 0
       uexact(i) = x(i)^2/2;
    elseif x(i) <= 0
       uexact(i) = -x(i)^2/2;
    end
end              

dt = 0.01*dx^2;       % Parameter of the iteration = stability condition
Nt = round(0.3/dt)    % Number of iterations of the Euler iteration                                 
u_new = zeros(Nx, 1);  
% Euler iterative process
for n = 1:Nt
   for i = 2:Nx-1
      % Define the function to be convolved   
      thetaeps(i-1)=theta1 * double(u(i) > 0) + theta2 * double(u(i) <0);
      % Upwind difference approximation for the spatial derivative
      ux(i) = max([-(u(i+1) - u(i))/dx,(u(i)-u(i-1))/dx,0]);
      % Central difference approximation for the 2nd order derivative 
      uxx(i) = (u(i+1)-2*u(i)+u(i-1))/(dx)^2;   
      % Update the solution using the finite difference scheme
     u_new(i) = u(i)-dt*(epsilon*u(i)+F(-uxx(i)))+dt*ff(i)/(power(epsilon+abs(ux(i))^2,thetaeps(i-1)/2));
    end
    % Boundary conditions 
    u_new(1) = -1/2;
    u_new(end) =1/2;
    % Update the solution
    u = u_new; 
    % Plot numerical solution at some iterations 
    figure(1)
     if mod(n, Nt/20) == 0
         plot(x, u,'LineWidth',1.5);
         hold on
         drawnow;
         pause(0.1); % Pause for less than half a second
     end
 end     
% Final plot
figure(2)
plot(x,v0(x)','LineWidth',3);
hold on
plot(x, u,'LineWidth',3);
plot(x, uexact,'LineWidth',3);
grid on
set(gca,'FontName','Times','FontSize',20);
xlabel('x','FontSize',20);
ylabel('u^n','FontSize',20);

%Compute the error
erro=norm(uexact-u',inf)