--- a/+scheme/Schrodinger1dCurve.m	Fri Feb 24 09:04:02 2017 +0100
+++ b/+scheme/Schrodinger1dCurve.m	Fri Feb 24 13:58:18 2017 +0100
@@ -10,6 +10,15 @@
         H % Discrete norm
         M % Derivative norm
         alpha
+        x_r
+        x_l
+        ddt_x_r
+        ddt_x_l
+        a
+        a_xi
+        Ji
+        t_up
+        x
         
         V_mat
         D1
@@ -24,16 +33,22 @@
     
     methods
         % Solving SE in the form u_t = i*u_xx +i*V on deforming 1D domain;
-        function obj = Schrodinger1dCurve(g,m,order,V,constJi)
+        function obj = Schrodinger1dCurve(g,order,boundaries,V,constJi)
             default_arg('V',0);
             default_arg('constJi',false)
             xilim={0 1};
+            m = N(g);
             if constJi
                 ops = sbp.D2Standard(m,xilim,order);
             else
                 ops = sbp.D4Variable(m,xilim,order);
             end
             
+            obj.x_l = boundaries{1};
+            obj.x_r = boundaries{2};
+            obj.ddt_x_l = boundaries{3};
+            obj.ddt_x_r = boundaries{4};
+            
             obj.xi=ops.x;
             obj.h=ops.h;
             obj.D2 = ops.D2;
@@ -47,15 +62,14 @@
             obj.d1_r = ops.d1_r;
             obj.grid = g;
             
-            
             if isa(V,'function_handle')
                 V_vec = V(obj.x);
             else
                 V_vec = obj.xi*0 + V;
             end
             
-            obj.V_mat = spdiags(V_vec,0,m,m);            
-            obj.D = @(a,a_xi,Ji) obj.d_fun(a, a_xi, Ji, constJi);           
+            obj.V_mat = spdiags(V_vec,0,m,m);
+            obj.D = @(t) obj.d_fun(t);
             obj.m = m;
             obj.order = order;
         end
@@ -69,11 +83,27 @@
         %       neighbour_scheme    is an instance of Scheme that should be interfaced to.
         %       neighbour_boundary  is a string specifying which boundary to interface to.
         
-        function [D] = d_fun(obj,a, a_xi , Ji , constJi)
-            if constJi
-                D= -0.5*(obj.D1*a - a_xi + a*obj.D1) + 1i*Ji*obj.D2 + 1i*obj.V_mat;
+        function [D] = d_fun(obj,t)
+   %         obj.update_vairables(t); In driscretization?
+            D= obj.Ji*(-0.5*(obj.D1*obj.a - obj.a_xi + obj.a*obj.D1) + 1i*obj.D2(diag(obj.Ji)) + 1i*obj.V_mat);
+            
+        end
+        
+        
+        function [] = variable_update(obj,t)
+            if (t == obj.t_up)
+                return
             else
-                D= -0.5*(obj.D1*a - a_xi + a*obj.D1) + 1i*obj.D2(diag(Ji)) + 1i*obj.V_mat;
+                x_r = obj.x_r(t);
+                x_l = obj.x_l(t);
+                ddt_x_r = obj. ddt_x_r(t);
+                ddt_x_l = obj.ddt_x_l(t);
+                obj.x = obj.xi*(x_r -x_l) + x_l;
+                obj.a = sparse(diag((-ddt_x_l*( x_r - x_l) - (obj.x-x_l)*(ddt_x_r-ddt_x_l))/(x_r-x_l)));            
+                
+                obj.Ji = sparse(diag(1./(x_r - x_l + 0*obj.x)));
+                obj.a_xi = sparse(diag(-1*(ddt_x_r - ddt_x_l + 0*obj.x)));
+                obj.t_up = t;
             end
         end
         
@@ -87,12 +117,12 @@
                 % Dirichlet boundary condition
                 case {'D','d','dirichlet'}
                     tau1 = s * 1i*d;
-                    tau2 = @(a) (-1*s*a(p,p) - abs(a(p,p)))/4*e;
-                    closure = @(a) obj.Hi*tau1*e' + obj.Hi*tau2(a)*e';
+                    tau2 = @(t) obj.Ji*(-1*s*obj.a(p,p) - abs(obj.a(p,p)))/4*e;
+                    closure = @(t) obj.Ji*obj.Hi*tau1*e' + obj.Hi*tau2(obj.a)*e';
                     
                     switch class(data)
                         case 'double'
-                            penalty = @(a) -(obj.Hi*tau1*data+obj.Hi*tau2(a)*data);
+                            penalty = @(t) -obj.Ji*(obj.Hi*tau1*data+obj.Hi*tau2(obj.a)*data);
                             %                      case 'function_handle'
                             %                           penalty = @(t)-obj.Hi*tau*data(t);
                         otherwise
@@ -108,18 +138,18 @@
         function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
             % u denotes the solution in the own domain
             % v denotes the solution in the neighbour domain
-                         [e_u,d_u,s_u,p_u] = obj.get_boundary_ops(boundary);
-                         [e_v,d_v,s_v,p_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
+            [e_u,d_u,s_u,p_u] = obj.get_boundary_ops(boundary);
+            [e_v,d_v,s_v,p_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
             
-                         a1 =  -s_u* 1/2 * 1i ;
-                         b1 =  a1';
-                         gamma = -1*s*a(p_u,p_u)/2;
+            a1 =  -s_u* 1/2 * 1i ;
+            b1 =  a1';
+            gamma = @(a) -s_u*a(p_u,p_u)/2*e_u;
             
-                         tau = b1*d_u;
-                         sig = -a1*e_u;
+            tau = b1*d_u;
+            sig = -a1*e_u;
             
-                         closure = @(a) obj.Hi * (tau*e_u' + sig*d_u' + gamma(a));
-                         penalty = @(a) obj.Hi * (-tau*e_v' - sig*d_v' - gamma(a));
+            closure = @(t) obj.Hi * (tau*e_u' + sig*d_u' + gamma(obj.a)*e_u');
+            penalty = @(t) obj.Hi * (-tau*e_v' - sig*d_v' - gamma(obj.a)*e_v');
         end
         
         % Ruturns the boundary ops and sign for the boundary specified by the string boundary.