Introduced function adapt_Δt, working for the case theta=0 (=keeping Δt constant).

If theta=!0, working but problem at with counter towards the end of the while loop.

eff_Drop.jl

@@ -3,6 +3,7 @@
 using Plots
 using DelimitedFiles
 using Interpolations
+using LinearAlgebra
 
 """
     analytical(x, x0, Nb, Tf, Di, Do, P)
@@ -153,6 +154,14 @@ function make_A(geometry,B,α_i, α_o, Δxi, Δxo, x, D_i, D_o, N, Nb, P, BEt)
     return A
 end
 
+function adapt_Δt(u, Δt, m_mat, Id, i, r, theta)
+    Heun=u[Id]-1/2*Δt*(m_mat[i]+m_mat[i+1])
+    Grad_desc=u[Id]-Δt*(m_mat[i])
+    δn=norm(Heun-Grad_desc)
+    Δt=Δt*(r/δn)^(theta/2)
+    return Δt
+end
+
 """
     solve_D_spherical(geometry,IC,AG,P,D_i,D_o,Δx,Δt,L,R,BP,x0,T)
 
@@ -277,9 +286,11 @@ function solve_D(geometry,IC,AG,B,P,D_i,D_o,Δx,Δt,L,R,BP,x0,T,perc,BE,tE)
     u_mat = zeros(Nt+1, length(u0))
     u_mat[1, :] = u0
     m_mat[1]=m(u0,Δxm)
+    # println(m_mat[1])
     u = A\b
     u_mat[2, :] = u
     m_mat[2]=m(u,Δxm)
+    # println(m_mat[2])
 
 
     #check for mass conservation
@@ -317,6 +328,22 @@ function solve_D(geometry,IC,AG,B,P,D_i,D_o,Δx,Δt,L,R,BP,x0,T,perc,BE,tE)
         u_mat[i, :] = u;
 
         m_mat[i]=m(u,Δxm)
+        # println(m_mat[i])
+
+        
+        #solution for the next time point if delta_t was kept constant
+        # and delta_t <- adapt_Δt
+        if t<T
+            b = calc_b(geometry,B, u, α_i(D_i,Δt,Δxi), α_o(D_o,Δt,Δxo), Δxi, Δxo,
+            x, Nb)
+            if B=="E"
+                A=make_A(geometry,B,α_i(D_i,Δt,Δxi), α_i(D_i,Δt,Δxo), Δxi, Δxo, x, D_i, D_o, length(x), Nb, P, BE_t(Δt*(i)))
+            end
+            u_next=A\b
+            m_mat[i+1]=m(u_next,Δxm)
+            # println(m_mat[i+1])
+            # Δt=adapt_Δt(u, Δt, m_mat, Id, i, 0.5, 0.001)
+        end
 
         t_vec[i]=t        
         #Mass conservation and recovery time for Neumann BC
@@ -335,11 +362,13 @@ function solve_D(geometry,IC,AG,B,P,D_i,D_o,Δx,Δt,L,R,BP,x0,T,perc,BE,tE)
                 t_perc = t
             end
         end
-        i+=1
         t+=Δt
         t=round(t, digits=length(string(Δt))-2)
+        # println(i)
+        # println(Δt)
+        # println(t)
+        i+=1
     end
 
     return A, u_mat, x, b_mat, Nb, t_perc, α_i(D_i,Δt,Δxi), α_o(D_o,Δt,Δxo), Mm, t_vec, Ni, No, m_mat
-end
-
+end
\ No newline at end of file