diff --git a/@bndfun/bndfun.m b/@bndfun/bndfun.m
index daa5b6cdff..af1bf666d4 100644
--- a/@bndfun/bndfun.m
+++ b/@bndfun/bndfun.m
@@ -77,7 +77,7 @@
                 pref = chebfunpref(pref);
             end
 
-            data = parseDataInputs(data, pref);
+            data = bndfun.parseDataInputs(data, pref);
 
             % Check the domain input.
             if ( ~all(size(data.domain) == [1, 2]) || (diff(data.domain) <= 0) )
@@ -120,12 +120,19 @@
         f = make(varargin);
         
     end
-    
-end
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% METHODS IMPLEMENTED IN THIS FILE:
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    %% PRIVATE STATIC METHODS:
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    methods ( Access = private, Static = true )
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Class-related functions: private utilities for this m-file.
+%% Note: temporarily (?) made private static methods for Octave
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
 
 function data = parseDataInputs(data, pref)
 %PARSEDATAINPUTS   Parse inputs from the DATA structure and assign defaults.
@@ -139,4 +146,7 @@
     data.hscale = norm(data.domain, inf);
 end
 
+end
+
+    end  % methods
 end
diff --git a/@chebfun/chebfun.m b/@chebfun/chebfun.m
index 3370f33108..a61b9caf0c 100644
--- a/@chebfun/chebfun.m
+++ b/@chebfun/chebfun.m
@@ -136,7 +136,7 @@
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     %% CLASS PROPERTIES:
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    properties (Access = public)
+    properties (Access = private)
         % DOMAIN of definition of a CHEBFUN object. If K = length(F.DOMAIN) is
         % greater than 1 then the CHEBFUN is referred to as a "piecewise".
         % CHEBFUN. The first and last values of this vector define the left and
@@ -144,8 +144,15 @@
         % locations of the interior breakpoints that define the domains of the
         % individual FUN objects comprising the CHEBFUN. The entries in this
         % vector should be strictly increasing.
-        domain              % (1x(K+1) double)
-
+        mydomain              % (1x(K+1) double)
+		
+		
+		% POINTVALUES Values of the function at the break points.
+        %pointValues = [];      % (1 x (K+1) double)
+		
+		
+    end
+    properties (Access = public)
         % FUNS is a cell array containing the FUN objects that comprise a
         % piecewise CHEBFUN. The kth entry in this cell is the FUN defining
         % the representation used by the CHEBFUN object on the open interval
@@ -163,7 +170,7 @@
         % (Mxinf) arrays. This difference is only behavioral; the other
         % properties described above are _NOT_ stored differently if this flag
         % is set.)
-        isTransposed = 0;   % (logical)
+        isTransposed = 0;   % (logical) ---!!!class field---
     end
     
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -200,7 +207,7 @@
             end
                        
             % Parse inputs:
-            [op, dom, data, pref, flags] = parseInputs(varargin{:});
+            [op, dom, data, pref, flags] = chebfun.parseInputs(varargin{:});
                         
             if ( flags.done )
                 % An update was performed. Exit gracefully:
@@ -216,14 +223,14 @@
                 % Construct from function_handle, numeric, or string input:
                 
                 % Call the main constructor:
-                [f.funs, f.domain] = chebfun.constructor(op, dom, data, pref);
+                [f.funs, f.mydomain] = chebfun.constructor(op, dom, data, pref);
                 
                 if ( flags.doubleLength )
                     % Using the length of f.funs{1} is okay because the
                     % 'doubleLength' flag is mutually exclusive with 'splitting
                     % on'.
                     pref.techPrefs.fixedLength = 2*length(f.funs{1}) - 1;
-                    [f.funs, f.domain] = chebfun.constructor(op, dom, data, pref);
+                    [f.funs, f.mydomain] = chebfun.constructor(op, dom, data, pref);
                 end
 
                 % Update values at breakpoints (first row of f.pointValues):
@@ -242,6 +249,41 @@
             end
 
         end
+		
+		
+		%{
+        # in case pointValues ​​is a private class field		
+        function set_private_point_values(f, array)
+			f.pointValues = array;		
+        end
+        
+
+        function f = set_private_point_values(f, varargin)
+
+          disp("set private property here!")
+          %approach 1!
+    
+          prop = varargin{1};
+          val = varargin{2};
+          while (numel (varargin) > 1)
+            prop = varargin{1};
+            val = varargin{2};
+            varargin(1:2) = [];
+    
+            f.pointValues = val(:); % force row vector
+          endwhile
+          disp("set result")
+          f.pointValues
+          disp("set result end")
+    
+        endfunction
+		
+		function res = get_private_point_values(f)
+            disp("get here!")
+			res = f.pointValues;		
+		end
+		
+        %}
         
     end
     
@@ -442,14 +484,14 @@
         % Parse inputs to PLOT. Extract 'lineWidth', etc.
         [lineStyle, pointStyle, jumpStyle, deltaStyle, out] = ...
             parsePlotStyle(varargin)
-    end
 
-end
 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Class-related functions: private utilities for this m-file.
+%% Note: temporarily (?) made private static methods for Octave
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
 function op = str2op(op)
     % Convert string inputs to either numeric format or function_handles.
     sop = str2num(op); %#ok<ST2NM> % STR2DOUBLE doesn't support str2double('pi')
@@ -779,10 +821,10 @@
     function op = parseOp(op)
         % Convert string input to function_handle:
         if ( ischar(op) )
-            op = str2op(op);
+            op = chebfun.str2op(op);
         end
         if ( doVectorCheck && isa(op, 'function_handle') )
-            op = vectorCheck(op, dom, vectorize);
+            op = chebfun.vectorCheck(op, dom, vectorize);
         end
         if ( isa(op, 'chebfun') )
             if ( op.isTransposed )
@@ -944,7 +986,6 @@
     
 end
 
-end
 
 function g = vec(op, y)
 %VEC  Vectorize a function or string expression.
@@ -984,3 +1025,10 @@
     end
 
 end
+
+end
+
+
+
+    end  % methods
+end
diff --git a/@chebfun/compose.m b/@chebfun/compose.m
index b02ea8ebfb..0b27c279fc 100644
--- a/@chebfun/compose.m
+++ b/@chebfun/compose.m
@@ -221,6 +221,7 @@
     
     if ( numel(f) == 1 )
         % Array-valued CHEBFUN case:
+        %disp(op)
         f = columnCompose(f, op, g, pref, opIsBinary);
     else
         % QUASIMATRIX case:
@@ -239,30 +240,41 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-function f = columnCompose(f, op, g, pref, opIsBinary)
+function f = columnCompose(obj, op, g, pref, opIsBinary)
+
+
 
 %% Initialise:
 
 % Initialise pointValues.
+f = chebfun();
+%disp("in columnCompose from compose.m!!!!!")
 
 if ( opIsBinary )
     % Call OVERLAP() if we are composing two CHEBFUN inputs with a binary op:
-    [f, g] = overlap(f, g);
-    pointValues = feval(op, f.pointValues, g.pointValues);
+    [f, g] = overlap(obj, g);
+    pointValues = feval(op, obj.pointValues, g.pointValues);
     newPointValues = pointValues(1,:);
 else
-    pointValues = feval(op, f.pointValues);
+    pointValues = feval(op, obj.pointValues);
     newPointValues = pointValues(1,:);
 end
 
 % Number of piecewise intervals in f:
-numInts = numel(f.domain) - 1;
+%disp(obj.mydomain)
+numInts = numel(obj.mydomain) - 1;
+%disp(numInts)
 
 % Initialise storage for the output FUN cell:
 newFuns = {};
 
-% Initialise new domain vector:
-newDom = f.domain(1);
+% Initialise new domain vector
+
+temp_dom = obj.mydomain;
+newDom = temp_dom(1);
+%newDom = get(f, "domain")(1); %%%% working 
+%newDom = f.mydomain(1); %%%%!!!!!! error: can't perform indexing operation on array of chebfun objects
+
 
 if ( pref.splitting) % Is splitting on?
     % Set the maximum length (i.e., number of sample points for CHEBTECH):
@@ -280,19 +292,22 @@
 % This line assumes that the compose method for f is the one that gets called
 % ultimately instead of the one for g (which is true currently but could
 % change).
-pref.tech = get(f.funs{1}, 'tech');
+temp_funs = obj.funs;
+pref.tech = get(temp_funs{1}, 'tech');
 
 % Suppress growing vector Mlint warnings (which are inevitable here):
 %#ok<*AGROW>
 
 %% Loop through each interval:
 for k = 1:numInts
+
+    %disp(k)
     
     % Attempt to generate FUN objects using FUN/COMPOSE().
     if ( isempty(g) )
-        newFun = compose(f.funs{k}, op, [], [], pref);
+        newFun = compose(temp_funs{k}, op, [], [], pref);
     else
-        newFun = compose(f.funs{k}, op, g.funs{k}, [], pref);
+        newFun = compose(temp_funs{k}, op, g.funs{k}, [], pref);
     end
     isHappy = get(newFun, 'ishappy');
 
@@ -314,7 +329,7 @@
         % Store new FUN in cell array:
         newFuns = [newFuns, {newFun}];                   
         % Store new ends:
-        newDom = [newDom, f.domain(k+1)];
+        newDom = [newDom, temp_dom(k+1)];
         % Store new pointValues: (Note, will only be a matrix - not a tensor)
         if ( isempty(g) )
             newPointValues = [newPointValues ; pointValues(k+1,:)];
@@ -325,7 +340,7 @@
     elseif ( pref.splitting )
 
         % If not happy and splitting is on, get a CHEBFUN for that subinterval:
-        domk = f.domain(k:k+1);
+        domk = temp_dom(k:k+1);
         if ( opIsBinary )
             newChebfun = chebfun(@(x) feval(op, feval(f, x), feval(g, x)), ...
                 domk, pref);
@@ -359,9 +374,11 @@
 %% Prepare output:
 
 % Put the FUN cell, domain, and pointValues back into a CHEBFUN:
-f.funs = newFuns;
-f.domain = newDom;
+
 f.pointValues = newPointValues;
+%f = set_private_point_values(f, "pointValues", newPointValues);
+f.funs = newFuns;
+f.mydomain = newDom;
 
 end
 
diff --git a/@chebfun/constructor.m b/@chebfun/constructor.m
index 4d2d62a4e6..e589fe03e6 100644
--- a/@chebfun/constructor.m
+++ b/@chebfun/constructor.m
@@ -124,6 +124,10 @@
     end
 end
 
+%disp("in constructor funs")
+%iscell(funs)
+%disp(funs)
+
 end
 
 %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%  SPLITTING ON %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/@chebfun/disp.m b/@chebfun/disp.m
index 4ede4cfb53..fe8595f9d1 100644
--- a/@chebfun/disp.m
+++ b/@chebfun/disp.m
@@ -7,7 +7,16 @@ function disp(f)
 % See http://www.chebfun.org/ for Chebfun information.
 
 % If the 'format loose' setting is enabled, we print additional linebreaks:
-loose = strcmp(get(0, 'FormatSpacing'), 'loose');
+if (exist('OCTAVE_VERSION', 'builtin') )
+    if (compare_versions(OCTAVE_VERSION(), '4.3.0', '>='))
+        [fmt, spacing] = format();
+        loose = strcmp(spacing, 'loose');
+    else
+        loose = eval('! __compactformat__ ()');
+    end
+else
+    loose = strcmp(get(0, 'FormatSpacing'), 'loose');
+end
 
 s = '';
 
@@ -70,9 +79,10 @@ function disp(f)
 s = [s, sprintf('\n       interval       length     endpoint values %s\n', extraItem)];
 len = zeros(numFuns, 1);
 for j = 1:numFuns
-    len(j) = length(f.funs{j});
+    ffuns = f.funs;
+    len(j) = length(ffuns{j});
 
-    if ( ~isreal(f.funs{j}) )
+    if ( ~isreal(ffuns{j}) )
         % For complex-valued funs, we don't display the values.
 
         % Print information to screen:
@@ -82,7 +92,7 @@ function disp(f)
     else
 
         % Grab values at endpoints:
-        endvals = [get(f.funs{j}, 'lval'), get(f.funs{j}, 'rval')];
+        endvals = [get(ffuns{j}, 'lval'), get(ffuns{j}, 'rval')];
 
         % Tweak the endpoint values: (This prevents -0 and +0)
         if ( ~any(isnan(endvals)) )
@@ -90,8 +100,9 @@ function disp(f)
         end
 
         % Print information to screen:
+        temp = f.domain;   % TODO: improve chaining in subsref?
         s = [s, sprintf('[%8.2g,%8.2g]   %6i  %8.2g %8.2g %s\n', ...
-            f.domain(j), f.domain(j+1), len(j), endvals, extraData{j})];
+            temp(j), temp(j+1), len(j), endvals, extraData{j})];
 
     end
 end
diff --git a/@chebfun/dispData.m b/@chebfun/dispData.m
index 0c15c7a529..bb7e5a8c5e 100644
--- a/@chebfun/dispData.m
+++ b/@chebfun/dispData.m
@@ -18,7 +18,8 @@
 
 % Loop over each FUN:
 for j = 1:numFuns
-    infoJ = dispData(f.funs{j});
+    temp = f.funs;
+    infoJ = dispData(temp{j});
     if ( ~isempty(infoJ) )
         numInfo = numel(infoJ);
         for k = 1:numInfo
diff --git a/@chebfun/display.m b/@chebfun/display.m
index 5533ec0561..f179d8f857 100644
--- a/@chebfun/display.m
+++ b/@chebfun/display.m
@@ -13,7 +13,18 @@ function display(X)
 % Copyright 2017 by The University of Oxford and The Chebfun Developers. 
 % See http://www.chebfun.org/ for Chebfun information.
 
-if ( isequal(get(0, 'FormatSpacing'), 'compact') )
+if (exist('OCTAVE_VERSION', 'builtin') )
+    if (compare_versions(OCTAVE_VERSION(), '4.3.0', '>='))
+        [fmt, spacing] = format();
+        loose = strcmp(spacing, 'loose');
+    else
+        loose = eval('! __compactformat__ ()');
+    end
+else
+    loose = strcmp(get(0, 'FormatSpacing'), 'loose');
+end
+
+if ( ~ loose )
     disp([inputname(1), ' =']);
     disp(X);
 else
diff --git a/@chebfun/domain.m b/@chebfun/domain.m
index 9afd3bb091..e8b3f6694d 100644
--- a/@chebfun/domain.m
+++ b/@chebfun/domain.m
@@ -11,29 +11,31 @@
 
 if ( nargout == 2 )
     % Return the end points as two scalars:
-    dom = f(1).domain([1 end]);
+    dom = f(1).mydomain([1 end]);
     A = dom(1);
     B = dom(2);
 
 elseif ( nargin == 2 )
     % Return the end points as a vector
     if ( strcmpi(flag, 'ends') )
-        A = f(1).domain([1 end]);
+        A = f(1).mydomain([1 end]);
     else
         error('CHEBFUN:CHEBFUN:domain:unknown', 'Unexpected input.');
     end
     
 elseif ( numel(f) == 1 )
     % CHEBFUN case:
-    A = f.domain;
+    A = f.mydomain;
     
 else
     % Merge the domains of columns in a quasimatrix:
     dom = cell(1, numel(f));
     for k = 1:numel(f)
-        dom{k} = f(k).domain;
+        dom{k} = f(k).mydomain;
     end
-    A = domain.merge(dom{:});
+    A = mydomain.merge(dom{:});
 end
 
 end
+
+
diff --git a/@chebfun/feval.m b/@chebfun/feval.m
index ad5cae3988..e2a43faf2f 100644
--- a/@chebfun/feval.m
+++ b/@chebfun/feval.m
@@ -201,4 +201,4 @@
         ones(1, numCols)).'), sizefx);
 end
 
-end
+end
\ No newline at end of file
diff --git a/@chebfun/get.m b/@chebfun/get.m
index 6dd5b40379..0744db7758 100644
--- a/@chebfun/get.m
+++ b/@chebfun/get.m
@@ -100,14 +100,6 @@
         out = hscale(f);
     case 'ishappy'
         out = ishappy(f);
-    case fieldnames(f)
-        out = cell(1, numel(f));
-        for k = 1:1:numel(f)
-            out{k} = f(k).(prop);
-        end
-        if ( numel(out) == 1 )
-            out = out{1};
-        end
     case 'lval'
         dom = domain(f);
         out = feval(f, dom(1));
@@ -167,6 +159,14 @@ case fieldnames(f)
         warning('CHEBFUN:CHEBFUN:get:imps', ...
             '''imps'' property is deprecated.  Use ''deltas'' instead.');
         out = get(f, 'deltas');
+    case fieldnames(f)
+        out = cell(1, numel(f));
+        for k = 1:1:numel(f)
+            out{k} = f(k).(prop);
+        end
+        if ( numel(out) == 1 )
+            out = out{1};
+        end
     otherwise
         error('CHEBFUN:CHEBFUN:get:badProp', ...
             '''%s'' is not a recognized CHEBFUN property name.', prop);
diff --git a/@chebfun/hscale.m b/@chebfun/hscale.m
index c8d2aee543..a91d3bc8f5 100644
--- a/@chebfun/hscale.m
+++ b/@chebfun/hscale.m
@@ -7,7 +7,7 @@
 % See http://www.chebfun.org/ for Chebfun information.
 
 % Compute INF norm of the domain:
-h = norm(f(1).domain([1, end]), inf);
+h = norm(f(1).mydomain([1, end]), inf);
 
 % Unbounded domains are defined to have hscale = 1:
 if ( isinf(h) )
diff --git a/@chebfun/isempty.m b/@chebfun/isempty.m
index b0dc573d74..8491399a23 100644
--- a/@chebfun/isempty.m
+++ b/@chebfun/isempty.m
@@ -7,7 +7,8 @@
 % See http://www.chebfun.org/ for Chebfun information.
 
 % If there is a .fun and f.funs(1) is not empty the f is not empty.
-if ( numel(f) > 1 || ( (numel(f.funs) > 0) && ~isempty(f.funs{1}) ) )
+temp = f.funs;
+if ( numel(f) > 1 || ( (numel(temp) > 0) && ~isempty(temp{1}) ) )
     out = false;
 else
     out = true;
diff --git a/@chebfun/mtimes.m b/@chebfun/mtimes.m
index 4704e2b9d5..48d4cf99be 100644
--- a/@chebfun/mtimes.m
+++ b/@chebfun/mtimes.m
@@ -88,7 +88,7 @@
 
 elseif ( isnumeric(f) && gIsChebfun )       % double * CHEBFUN
 
-        f = mtimes(g.', f.').';
+        f = mtimes(g, f);
 
 elseif ( fIsChebfun )                       % CHEBFUN * ??? 
 
@@ -102,4 +102,4 @@
 
 end
 
-end
+end
\ No newline at end of file
diff --git a/@chebfun/numColumns.m b/@chebfun/numColumns.m
index 0e89a49df9..fe772dd254 100644
--- a/@chebfun/numColumns.m
+++ b/@chebfun/numColumns.m
@@ -13,7 +13,8 @@
     out = 0;
 elseif ( numel(f) == 1 )
     % Possible array-valued:
-    out = size(f.funs{1}, 2);
+    temp = f.funs;
+    out = size(temp{1}, 2);
 else
     % A quasimatrix:
     out = numel(f);
diff --git a/@chebfun/plus.m b/@chebfun/plus.m
index 51cd01534a..7699fa4c35 100644
--- a/@chebfun/plus.m
+++ b/@chebfun/plus.m
@@ -1,4 +1,4 @@
-function f = plus(f, g)
+function f = plus(obj, g)
 %+   CHEBFUN plus.
 %   F + G adds CHEBFUNs F and G, or a scalar to a CHEBFUN if either F or G is a
 %   scalar.
@@ -10,11 +10,12 @@
 
 % Copyright 2017 by The University of Oxford and The Chebfun Developers.
 % See http://www.chebfun.org/ for Chebfun information.
+f = chebfun();
 
-if ( ~isa(f, 'chebfun') )   % ??? + CHEBFUN
+if ( ~isa(obj, 'chebfun') )   % ??? + CHEBFUN
     
     % Ensure CHEBFUN is the first input:
-    f = plus(g, f);
+    f = plus(g, obj);
     return
     
 elseif ( isempty(g) )       % CHEBFUN + []
@@ -23,31 +24,43 @@
     return
     
 elseif ( isnumeric(g) )     % CHEBFUN + double
+
+    temp_funs = obj.funs;
+    temp_pointValues = obj.pointValues;
+    temp_domain = obj.mydomain;
     
-    if ( numel(f) == 1 )
+    if ( numel(obj) == 1 )
         % Array-valued case:
         
         % Transpose g if f is transposed:
-        if ( f(1).isTransposed )
+        if ( obj(1).isTransposed )
             g = g.';
         end
         
         % Add g to the FUNs:
-        for k = 1:numel(f.funs)
-            f.funs{k} = f.funs{k} + g;
+	
+        for k = 1:numel(obj.funs)
+            f.funs{k} = temp_funs{k} + g;
         end
+	% todo: reassign?
         % Add g to the pointValues:
-        if ( (size(f.pointValues, 2) == 1) &&  (numColumns(f) > 1) )
-            f.pointValues = repmat(f.pointValues, 1, size(g, 2)); % Allow expansion in f.
+        if ( (size(temp_pointValues, 2) == 1) &&  (numColumns(f) > 1) )
+            f.pointValues = repmat(temp_pointValues, 1, size(g, 2)); % Allow expansion in f.
         end
         if ( size(g, 2) > 1 )
-            g = repmat(g, length(f.domain), 1);             % Allow expansion in g.
+            g = repmat(g, length(obj.domain), 1);             % Allow expansion in g.
         end
-        f.pointValues = f.pointValues + g;
+        %disp("before set")
+        %disp(obj.pointValues)
+        
+        f.pointValues = temp_pointValues + g; 
+        %f = set_private_point_values(f, "pointValues", temp_pointValues + g); 
+        %disp("after set")
+        %f.pointValues = temp_pointValues + g;
     else
         % Quasimatrix case:
         
-        numCols = numel(f);
+        numCols = numel(obj);
         % Promote g if required:
         if ( isscalar(g) )
             g = repmat(g, 1, numCols);
@@ -71,61 +84,65 @@
         ['Undefined function ''plus'' for input arguments of type %s ' ...
         'and %s.'], class(f), class(g));
     
-elseif ( isempty(f) )         % empty CHEBFUN + CHEBFUN
+elseif ( isempty(obj) )         % empty CHEBFUN + CHEBFUN
     
     % Nothing to do. (Return empty CHEBFUN as output).
     return
     
 else                          % CHEBFUN + CHEBFUN
     
+    temp_funs = obj.funs;
+    temp_pointValues = obj.pointValues;
+    temp_domain = obj.mydomain;
+    
     % Check to see if one CHEBFUN is transposed:
     if ( xor(f(1).isTransposed, g(1).isTransposed) )
         error('CHEBFUN:CHEBFUN:plus:matdim', ...
             'Matrix dimensions must agree. (One input is transposed).');
     end
         
-    dimCheck(f, g);
+    dimCheck(obj, g);
     
-    if ( numel(f) == 1 && numel(g) == 1 )
+    if ( numel(obj) == 1 && numel(g) == 1 )
         % CHEBFUN case:
         
         % If one of the two CHEBFUNs uses a PERIODICTECH reprensetation, 
         % cast it to a NONPERIODICTECH.
-        if ( ~isPeriodicTech(f.funs{1}) && isPeriodicTech(g.funs{1}) )
+        if ( ~isPeriodicTech(temp_funs{1}) && isPeriodicTech(g.funs{1}) )
             g = chebfun(g, g.domain, 'tech', get(f.funs{1}, 'tech'));
-        elseif ( isPeriodicTech(f.funs{1}) && ~isPeriodicTech(g.funs{1}) )
-            f = chebfun(f, f.domain, 'tech', get(g.funs{1}, 'tech'));
+        elseif ( isPeriodicTech(temp_funs{1}) && ~isPeriodicTech(g.funs{1}) )
+            f = chebfun(obj, obj.domain, 'tech', get(g.funs{1}, 'tech'));
         end
         
         % Overlap the CHEBFUN objects:
-        [f, g] = overlap(f, g);
+        [f, g] = overlap(obj, g);
         
         % Add the pointValues:
-        f.pointValues = f.pointValues + g.pointValues;
+        f.pointValues = temp_pointValues + g.pointValues;
         % Add the FUNs:
-        for k = 1:numel(f.funs)
-            f.funs{k} = f.funs{k} + g.funs{k};
+        for k = 1:numel(temp_funs)
+            f.funs{k} = temp_funs{k} + g.funs{k};
         end
 
     else
         % QUASIMATRIX case:
         
         % Convert to cell for simplicity
-        f = cheb2cell(f);
+        f = cheb2cell(obj);
         g = cheb2cell(g);
         
         % Loop over the columns:
-        if ( numel(f) == 1 )
+        if ( numel(obj) == 1 )
             for k = numel(g):-1:1
-                h(k) = f{1} + g{k};
+                h(k) = obj{1} + g{k};
             end
         elseif ( numel(g) == 1 )
             for k = numel(f):-1:1
-                h(k) = f{k} + g{1};
+                h(k) = obj{k} + g{1};
             end
         else % numel(f) = numel(g)
-            for k = numel(f):-1:1
-                h(k) = f{k} + g{k};
+            for k = numel(obj):-1:1
+                h(k) = obj{k} + g{k};
             end
         end
         f = h;
@@ -134,7 +151,14 @@
 
 end
 
+f.mydomain = temp_domain;
+
+%disp(f);
+%disp(f.funs);
+%disp(f.mydomain);
+%disp(f.pointValues);
+
 % Set small breakpoint values to zero:
 f = thresholdBreakpointValues(f);
 
-end
+end
\ No newline at end of file
diff --git a/@chebfun/subsasgn.m b/@chebfun/subsasgn.m
index 972e92226a..75b4872d90 100644
--- a/@chebfun/subsasgn.m
+++ b/@chebfun/subsasgn.m
@@ -25,6 +25,13 @@
 
 idx = index(1).subs;
 
+% debugging:
+disp("we're in @chebfun/subsasgn...")
+index
+idx
+val
+
+
 switch index(1).type
     case '.'
 
@@ -54,4 +61,4 @@
             ['??? Unexpected index.type of ' index(1).type]);
 end
 
-end
+end
\ No newline at end of file
diff --git a/@chebfun/subsref.m b/@chebfun/subsref.m
index 90c9154f44..9bd2aa1fb8 100644
--- a/@chebfun/subsref.m
+++ b/@chebfun/subsref.m
@@ -38,6 +38,7 @@
 
 %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%% FEVAL / COMPOSE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     case '()'
+        %disp('()');
         
         % Deal with row CHEBFUN objects:
         isTransposed = f(1).isTransposed;
@@ -144,6 +145,7 @@
     
 %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% GET %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     case '.'
+        %disp('.');
 
         % Call GET() for .PROP access.
         out = get(f, idx);
@@ -155,6 +157,7 @@
 
 %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% RESTRICT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     case '{}'
+        %disp('{}');
 
         if ( length(idx) == 1 )
             if ( isequal(idx{1}, ':') )
diff --git a/@chebfun/times.m b/@chebfun/times.m
index 7c45102e6e..538c102eec 100644
--- a/@chebfun/times.m
+++ b/@chebfun/times.m
@@ -125,4 +125,4 @@
 % Set small breakpoint values to zero:
 f = thresholdBreakpointValues(f);
 
-end
+end
\ No newline at end of file
diff --git a/@classicfun/classicfun.m b/@classicfun/classicfun.m
index 939442c2e0..f7e84d15bd 100644
--- a/@classicfun/classicfun.m
+++ b/@classicfun/classicfun.m
@@ -92,7 +92,7 @@
                 pref = chebfunpref(pref);
             end
 
-            data = parseDataInputs(data, pref);
+            data = classicfun.parseDataInputs(data, pref);
 
             % Call constructor depending on domain:
             if ( ~any(isinf(data.domain)) )
@@ -164,12 +164,18 @@
         end
 
     end
-end
+
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% METHODS IMPLEMENTED IN THIS FILE:
+%% Class-related functions: private utilities for this m-file.
+%% Note: temporarily (?) made private static methods for Octave
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    %% PRIVATE STATIC METHODS:
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    methods ( Access = private, Static = true )
+
 function data = parseDataInputs(data, pref)
 %PARSEDATAINPUTS   Parse inputs from the DATA structure and assign defaults.
 
@@ -177,4 +183,8 @@
     data.domain = pref.domain;
 end
 
+end
+
+
+    end  % methods
 end
diff --git a/@domain/cumsum.m b/@domain/cumsum.m
deleted file mode 100644
index 3669f034ab..0000000000
--- a/@domain/cumsum.m
+++ /dev/null
@@ -1,10 +0,0 @@
-function J = cumsum(d, varargin)
-%CUMSUM    Indefinite integration operator.
-%   This function is deprecated. Use OPERATORBLOCK.CUMSUM.
-
-% Copyright 2017 by The University of Oxford and The Chebfun Developers. 
-% See http://www.chebfun.org/ for Chebfun information.
-
-J = linop( operatorBlock.cumsum(double(d), varargin{:}) );
-
-end
diff --git a/@domain/diag.m b/@domain/diag.m
deleted file mode 100644
index ddfb24e5ec..0000000000
--- a/@domain/diag.m
+++ /dev/null
@@ -1,10 +0,0 @@
-function F = diag(f, d)
-%DIAG      Pointwise multiplication operator.
-%   This function is deprecated. Use OPERATORBLOCK.MULT.
-
-% Copyright 2017 by The University of Oxford and The Chebfun Developers.
-% See http://www.chebfun.org/ for Chebfun information.
-
-F = linop( operatorBlock.mult(f, double(d)) );
-
-end
diff --git a/@domain/diff.m b/@domain/diff.m
deleted file mode 100644
index be2c4e2aee..0000000000
--- a/@domain/diff.m
+++ /dev/null
@@ -1,10 +0,0 @@
-function F = diff(d, varargin)
-%DIFF      Differentiation operator.
-%   This function is deprecated. Use OPERATORBLOCK.DIFF.
-
-% Copyright 2017 by The University of Oxford and The Chebfun Developers.
-% See http://www.chebfun.org/ for Chebfun information.
-
-F = linop( operatorBlock.diff(double(d), varargin{:}) );
-
-end
diff --git a/@domain/domain.m b/@domain/domain.m
deleted file mode 100644
index cc845a025c..0000000000
--- a/@domain/domain.m
+++ /dev/null
@@ -1,107 +0,0 @@
-classdef (InferiorClasses = {?chebfun}) domain < double
-%DOMAIN   Utility class for CHEBFUN. Mostly for backward compatibility.
-
-% Copyright 2017 by The University of Oxford and The Chebfun Developers.
-% See http://www.chebfun.org/ for Chebfun information.
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% DOMAIN Class Description:
-%
-% DOMAIN inherits from a standard Matlab DOUBLE. A domain object only
-% contains vector for the endpoints and breakpoints of the interval it
-% represents. This class is lightly documented, since it is mostly intended
-% for backward compatibility.
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    
-    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    %% CLASS PROPERTIES:
-    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    properties ( Access = public )
-        
-    end
-    
-    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    %% CLASS CONSTRUCTOR:
-    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    methods ( Access = public, Static = false )
-
-        function obj = domain(varargin)
-            %Constructor for the DOMIAN class.
-            
-            % Return an empty DOMAIN on null input:
-            if ( nargin == 0 )
-                data = [];
-            else
-                data = horzcat(varargin{:});
-            end
-            
-            % Create the domain:
-            obj = obj@double(data);
-            
-        end
-        
-    end     
-    
-    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    %% CLASS METHODS:
-    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    methods ( Access = public, Static = false )
-        
-        function out = validate(dom)
-            out = true;
-            if ( ~isnumeric(dom) || any(isnan(dom)) )
-                out = false;
-            elseif ( isempty(dom) )
-                return
-            elseif ( size(dom, 1) > 1 || size(dom, 2) < 2 )
-                out = false;
-            elseif ( any(diff(dom) <= 0) )
-                out = false;
-            end
-            if ( out == false && nargout == 0 )
-                error('CHEBFUN:DOMAIN:domain:invalid', ...
-                    'Ends must be a 1xM vector of ordered doubles.');
-            end
-        end
-        
-        function display(dom)
-            disp(double(dom))
-        end
-        
-        function out = subsref(d, s)
-            out = subsref@double(d, s);
-            out = double(out);
-        end
-        
-        function varargout = sprintf(varargin)
-            % This is required as built-in subsref does not know what to do with
-            % a DOMAIN object.
-            
-            varargin = domain2double(varargin{:});
-            
-            % Call built-in SPRINTF:
-            varargout{1:nargout} = sprintf(varargin{:});
-        end
-                   
-    end
-    
-    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    %% STATIC METHODS:
-    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    methods ( Access = public, Static = true )
-        
-        % Merge two domains:
-        newDom = merge(varargin)
-        
-        function varargin = toDouble(varargin)
-            % Cast DOMAIN to DOUBLE:
-            for k = 1:nargin
-                if ( isa(varargin{k}, 'domain') )
-                    varargin{k} = double(varargin{k});
-                end
-            end
-        end
-        
-    end
-    
-end
diff --git a/@domain/eye.m b/@domain/eye.m
deleted file mode 100644
index 8d21674a31..0000000000
--- a/@domain/eye.m
+++ /dev/null
@@ -1,10 +0,0 @@
-function I = eye(d)
-%EYE       Identity operator.
-%   This function is deprecated. Use OPERATORBLOCK.EYE.
-
-% Copyright 2017 by The University of Oxford and The Chebfun Developers.
-% See http://www.chebfun.org/ for Chebfun information.
-
-I = linop( operatorBlock.eye(double(d)) );
-
-end
diff --git a/@domain/feval.m b/@domain/feval.m
deleted file mode 100644
index 71a74868e1..0000000000
--- a/@domain/feval.m
+++ /dev/null
@@ -1,10 +0,0 @@
-function F = feval(d, loc, varargin)
-%FEVAL     Evaluation functional.
-%   This function is deprecated. Use FUNCTIONALBLOCK.FEVAL.
-
-% Copyright 2017 by The University of Oxford and The Chebfun Developers.
-% See http://www.chebfun.org/ for Chebfun information.
-
-F = linop( functionalBlock.feval(loc, double(d), varargin{:}) );
-
-end
diff --git a/@domain/fred.m b/@domain/fred.m
deleted file mode 100644
index c11cfc9ca2..0000000000
--- a/@domain/fred.m
+++ /dev/null
@@ -1,10 +0,0 @@
-function F = fred(k, d, varargin)
-%FRED      Fredholm integral operator.
-%   This function is deprecated. Use OPERATORBLOCK.FRED.
-
-% Copyright 2017 by The University of Oxford and The Chebfun Developers.
-% See http://www.chebfun.org/ for Chebfun information.
-
-F = linop( operatorBlock.fred(k, double(d), varargin{:}) );
-
-end
diff --git a/@domain/interp1.m b/@domain/interp1.m
deleted file mode 100644
index 36514a0f4b..0000000000
--- a/@domain/interp1.m
+++ /dev/null
@@ -1,45 +0,0 @@
-function p = interp1(x, y, method, dom)
-%INTERP1   CHEBFUN polynomial interpolant at any distribution of points.
-%   P = INTERP1(X, Y, D), where X and Y are vectors, returns the CHEBFUN P
-%   defined on the domain D([1, end]) corresponding to the polynomial
-%   interpolant through the data Y(j) at points X(j).
-%
-%   If Y is a matrix with more than one column then then Y(:,j) is taken as the
-%   value to be matched at X(j) and P is an array-valued CHEBFUN with each
-%   column corresponding to the appropriate interpolant.
-%
-%   EXAMPLE: The following commands plot the interpolant in 11 equispaced points
-%   on [-1, 1] through the famous Runge function:
-%       d = domain(-1, 1);
-%       ff = @(x) 1./(1+25*x.^2);
-%       x = linspace(d(1), d(2), 11);
-%       p = interp1(x, ff(x), d)
-%       f = chebfun(ff, d);
-%       plot(f, 'k', p, 'r-'), hold on, plot(x, ff(x), 'r.'), grid on
-%
-%   P = CHEBFUN.INTERP1(X, Y, METHOD, D) specifies alternate interpolation
-%   methods. The default is as described above. (Use an empty matrix [] to
-%   specify the default.) Available methods are:
-%       'linear'   - linear interpolation
-%       'spline'   - piecewise cubic spline interpolation (SPLINE)
-%       'pchip'    - shape-preserving piecewise cubic interpolation
-%       'cubic'    - same as 'pchip'
-%       'poly'     - polynomial interpolation, as described above
-%
-% See also SPLINE, PCHIP.
-
-% Copyright 2017 by The University of Oxford and The Chebfun Developers.
-% See http://www.chebfun.org/ for Chebfun information.
-
-if ( nargin == 3 )
-    dom = method;
-    method = [];
-end
-
-if ( ~isnumeric(y) )
-    y = feval(y, x);
-end
-
-p = chebfun.interp1(x, y, method, double(dom));
-
-end
diff --git a/@domain/merge.m b/@domain/merge.m
deleted file mode 100644
index f822543e77..0000000000
--- a/@domain/merge.m
+++ /dev/null
@@ -1,69 +0,0 @@
-function newDom = merge(varargin)
-%MERGE    Merge breakpoints (with a tolerance).
-%   MERGE(DOM1, DOM2, ..., DOMN, TOL) merges the breakpoints of the domains
-%   DOM1, ..., DOMN to a tolerance TOL, where DOM1, ..., DOMN are sorted vectors
-%   of real numbers. If the domains are not compatible, i.e.. the first and
-%   final entry of each DOM differ by more than TOL, then an error is returned.
-%
-%   MERGE(DOM1, DOM2, ..., DOMN) uses a tolerance of the largest magnitude entry
-%   in DOM1, ..., DOMN scaled by 10*eps.
-%
-%   MERGE() always returns an array of doubles, even if the inputs contain
-%   DOMAIN objects.
-%
-% See also WHICHDOMAIN, TWEAKDOMAIN, DOMAINCHECK.
-
-% Copyright 2017 by The University of Oxford and The Chebfun Developers.
-% See http://www.chebfun.org/ for Chebfun information.
-
-% Relabel the input variable:
-doms = varargin;
-% Make sure we have doubles, not domains:
-doms = cellfun(@(dom) double(dom), doms, 'UniformOutput', false);
-
-% Ignore empties:
-doms(cellfun(@isempty, doms)) = [];
-
-% Initialise newDom:
-newDom = [];
-
-% Choose a tolerance:
-if ( isscalar(varargin{end}) )
-    % One is given:
-    tol = varargin{end};
-    doms(end) = [];
-else
-    % 100*eps*hscale:
-    hscales = cellfun(@(d) norm(d, inf), doms);
-    hscales(isinf(hscales)) = 1;
-    tol = 100*eps*max(hscales);
-end
-% Check to see if the domains are compatible:
-ends = cellfun(@(dom) dom([1 end]), doms, 'UniformOutput', false);
-diffEnds = cell2mat(ends.') - repmat(ends{1}, numel(doms), 1);
-if ( any(diffEnds(:) > tol) )
-    error('CHEBFUN:DOMAIN:merge:incompat', 'Incompatible domains.');
-end
-
-j = 1;
-while ( ~isempty(doms) )
-    % Find the minimum remaining domain entry:
-    newDom(j) = min(cellfun(@min, doms)); %#ok<AGROW>
-    for k = numel(doms):-1:1
-        % Remove -inf
-        doms{k}(isinf(doms{k}) & doms{k} < 0) = [];
-        % Remove subsequent domains less than tol away:
-        doms{k}(doms{k} < newDom(j) + tol) = [];
-        if ( isempty(doms{k}) )
-            doms(k) = []; % Discard if this domain is now empty.
-        end
-        if ( isinf(newDom(j)) && newDom(j) > 0 )
-            % Stop if we get to +inf
-            doms = [];
-            break
-        end
-    end
-    j = j + 1;
-end 
-
-end
diff --git a/@domain/ode113.m b/@domain/ode113.m
deleted file mode 100644
index 0a9df398ee..0000000000
--- a/@domain/ode113.m
+++ /dev/null
@@ -1,17 +0,0 @@
-function varargout = ode113(varargin)
-%ODE113   Solve stiff differential equations and DAEs. Output a CHEBFUN.
-%   
-% This syntax is deprecated. Please use chebfun.ode113(...) instead.
-
-% Copyright 2017 by The University of Oxford and The Chebfun Developers. 
-% See http://www.chebfun.org/ for Chebfun information.
-
-warning('CHEBFUN:DOMAIN:ode113:deprecated', ...
-    ['Usage of ODE113 via the @DOMAIN class is deprecated and may be ', ...
-    'removed from future releases. Please use chebfun.ode113(...) instead.']);
-
-varargin = domain.toDouble(varargin{:});
-
-varargout{1:nargout} = chebfun.ode113(varargin{:});
-
-end
diff --git a/@domain/ode15s.m b/@domain/ode15s.m
deleted file mode 100644
index 8421ad9970..0000000000
--- a/@domain/ode15s.m
+++ /dev/null
@@ -1,17 +0,0 @@
-function varargout = ode15s(varargin)
-%ODE15s   Solve stiff differential equations and DAEs. Output a CHEBFUN.
-%   
-% This syntax is deprecated. Please use chebfun.ode15s(...) instead.
-
-% Copyright 2017 by The University of Oxford and The Chebfun Developers. 
-% See http://www.chebfun.org/ for Chebfun information.
-
-warning('CHEBFUN:DOMAIN:ode15s:deprecated', ...
-    ['Usage of ODE15S via the @DOMAIN class is deprecated and may be ', ...
-    'removed from future releases. Please use chebfun.ode15s(...) instead.']);
-
-varargin = domain.toDouble(varargin{:});
-
-varargout{1:nargout} = chebfun.ode15s(varargin{:});
-
-end
diff --git a/@domain/ode45.m b/@domain/ode45.m
deleted file mode 100644
index e4cbfb08da..0000000000
--- a/@domain/ode45.m
+++ /dev/null
@@ -1,17 +0,0 @@
-function varargout = ode45(varargin)
-%ODE45   Solve non-stiff ODEs. Output a CHEBFUN.
-%   
-% This syntax is deprecated. Please use chebfun.ode45(...) instead.
-
-% Copyright 2017 by The University of Oxford and The Chebfun Developers. 
-% See http://www.chebfun.org/ for Chebfun information.
-
-warning('CHEBFUN:DOMAIN:ode45:deprecated', ...
-    ['Usage of ODE45 via the @DOMAIN class is deprecated and may be ', ...
-    'removed from future releases. Please use chebfun.ode45(...) instead.']);
-
-varargin = domain.toDouble(varargin{:});
-
-varargout{1:nargout} = chebfun.ode45(varargin{:});
-
-end
diff --git a/@domain/poly.m b/@domain/poly.m
deleted file mode 100644
index 2571230a9f..0000000000
--- a/@domain/poly.m
+++ /dev/null
@@ -1,65 +0,0 @@
-function f = poly(v, d)
-%POLY   Convert roots to CHEBFUN.
-%   POLY(V, D), when V is a column vector and D is a domain, returns a CHEBFUN
-%   of degree length(V) and domain D whose roots are the elements of V.
-%
-% Example (a nonstandard way to construct chebpoly(100)):
-%   T100 = poly(chebpts(100,1),domain(-1,1)); T100 = T100/T100(1);
-%
-% See also CHEBFUN/POLY.
-
-% Copyright 2017 by The University of Oxford and The Chebfun Developers.
-% See http://www.chebfun.org/ for Chebfun information.
-
-% Deal with array input:
-if ( min(size(v)) > 1 )
-    f = cell(1, size(v, 2));
-    for k = 1:size(v, 2)
-        f{k} = poly(v(:,k), d);
-    end
-    f = horzcat(f{:});
-    return
-end
-
-d = double(d);
-N = length(v);
-
-% Remove infs, and return NaN if NaN present:
-v = v(~isinf(v));
-if ( any(isnan(v)) )
-    f = chebfun(NaN, d);
-    return
-end
-
-% Return empty CHEBFUN if v is empty:
-if ( N == 0 )
-    f = chebfun();
-    return
-end
-
-% Leja ordering:
-[ignored, j] = max(abs(v));
-z = v(j);
-v(j) = [];
-for k = 1:(N-1)
-    P = zeros(N - k, 1);
-    for l = 1:(N-k)
-        P(l) = prod(z - v(l));
-    end
-    [ignored, j] = max(abs(P));
-    z(k+1) = v(j);
-    v(j) = [];
-end
-v = z;
-
-% Evaluate at Chebyshev points:
-x = chebpts(N+1, d);
-p = ones(N+1, 1);
-for k = 1:N
-    p = p.*(x - v(k));
-end
-
-% Contruct the CHEBFUN:
-f = chebfun(p, d, 'tech', @chebtech2);
-
-end
diff --git a/@domain/polyfit.m b/@domain/polyfit.m
deleted file mode 100644
index 7165e51247..0000000000
--- a/@domain/polyfit.m
+++ /dev/null
@@ -1,52 +0,0 @@
-function f = polyfit(x, y, n, d)
-%POLYFIT   Polyfit discrete data and return a CHEBFUN object.
-%   F = POLYFIT(X, Y, N, D), where D is a DOMAIN object, returns a CHEBFUN F on
-%   the domain D([1, end]) which corresponds to the polynomial of degree N that
-%   fits the data (X, Y) in the least-squares sense. X should be a real-valued
-%   column vector and Y should be a matrix with size(Y,1) = size(X,1).
-%
-%   Note DOMAIN/POLYFIT does not not support more than one output argument in
-%   the way that MATLAB/POLYFIT does.
-%
-% See also CHEBFUN/POLYFIT.
-
-% Copyright 2017 by The University of Oxford and The Chebfun Developers.
-% See http://www.chebfun.org/ for Chebfun information.
-
-% Convert domain to a double:
-d = double(d);
-d = d([1, end]);
-
-% Align dimensions:
-if ( size(x, 2) > 1 )
-    if ( size(x, 1) > 1 )
-        error('CHEBFUN:DOMAIN:polyfit:xIn', ...
-            'X should be a real-valued column vector');
-    end
-    x = x.';
-end
-if ( size(y, 1) ~= size(x, 1) )
-    if ( size(y, 2) ~= size(x, 1) )
-        error('CHEBFUN:DOMAIN:polyfit:xIn', ...
-            'X and Y vectors must be the same size.');
-    end
-    y = y.';
-end
-
-% Make Chebyshev-Vandermonde matrix:
-% The code below is a faster version of 
-%       T = chebpoly(0:n, d); Tx = feval(T, x);
-m = numel(x)-1;
-Tx = zeros( m+1, n+1); 
-Tx(:,1) = ones(m+1,1);
-x_map = 2*(x-d(1))/(d(2)-d(1)) - 1;
-Tx(:,2) = x_map;
-for k = 2:n
-    Tx(:,k+1) = 2*x_map.*Tx(:,k) - Tx(:,k-1);
-end
-% Solve for coefficients (least squares)
-c = Tx\y;
-% Construct Chebfun:
-f = chebfun(c, d, 'coeffs');
-
-end
diff --git a/@domain/sum.m b/@domain/sum.m
deleted file mode 100644
index 2adcba6ef6..0000000000
--- a/@domain/sum.m
+++ /dev/null
@@ -1,10 +0,0 @@
-function S = sum(d)
-%SUM       Integration functional.
-%   This function is deprecated. Use FUNCTIONALBLOCK.SUM.
-
-% Copyright 2017 by The University of Oxford and The Chebfun Developers.
-% See http://www.chebfun.org/ for Chebfun information.
-
-S = linop( functionalBlock.sum(double(d)) );
-
-end
diff --git a/@domain/volt.m b/@domain/volt.m
deleted file mode 100644
index 20d8cc880d..0000000000
--- a/@domain/volt.m
+++ /dev/null
@@ -1,10 +0,0 @@
-function V = volt(k, d, varargin)
-%VOLT      Volterra integral operator.
-%   This function is deprecated. Use OPERATORBLOCK.VOLT.
-
-% Copyright 2017 by The University of Oxford and The Chebfun Developers.
-% See http://www.chebfun.org/ for Chebfun information.
-
-V = linop( operatorBlock.volt(k, double(d), varargin{:}) );
-
-end
diff --git a/@domain/zeros.m b/@domain/zeros.m
deleted file mode 100644
index 2f42e19e23..0000000000
--- a/@domain/zeros.m
+++ /dev/null
@@ -1,10 +0,0 @@
-function Z = zeros(d)
-%ZEROS     Zero operator.
-%   This function is deprecated. Use OPERATORBLOCK.ZEROS.
-
-% Copyright 2017 by The University of Oxford and The Chebfun Developers.
-% See http://www.chebfun.org/ for Chebfun information.
-
-Z = linop( operatorBlock.zeros(double(d)) );
-
-end
diff --git a/@smoothfun/smoothfun.m b/@smoothfun/smoothfun.m
index f2b0d4a974..869ce1839d 100644
--- a/@smoothfun/smoothfun.m
+++ b/@smoothfun/smoothfun.m
@@ -77,11 +77,17 @@
         
     end
     
-end
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    %% STATIC METHODS:
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    methods ( Access = public, Static = true )
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Class-related functions: private utilities for this m-file.
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Note: temporarily (?) made private static methods for Octave
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 function f = funqui(vals)
 %FUNQUI   Rational interpolant of equispaced data.
@@ -196,4 +202,9 @@
     end
 end
 
+end
+
+
+    end  % methods
+
 end
diff --git a/chebfun_MWE.m b/chebfun_MWE.m
new file mode 100644
index 0000000000..b1b5bb4eb2
--- /dev/null
+++ b/chebfun_MWE.m
@@ -0,0 +1,52 @@
+% A container for Octave Chebfun examples
+% As we try to port Chebfun to Octave, let's keep track of what works.
+
+
+disp("chebfun initialization: ");
+x = chebfun('x', [-10 10])
+
+disp("\nexamples of product : ");
+disp("2*x");
+disp(2*x);
+disp("\n2.*x");
+disp(2.*x);
+disp("\nx*2");
+disp(x*2);
+disp("\nx.*2");
+disp(x.*2);
+
+
+disp("\nan example of division : ");
+disp("x/2");
+disp(x/2);
+
+disp("\nan example of addition : ");
+disp("x + x");
+disp(x + x);
+
+disp("\ntranscendental functions examples : ");
+disp("sin(x)");
+disp(sin(x));
+disp("\ncos(x)");
+disp(cos(x));
+
+
+disp("\nanother more complicated examples : ");
+disp("x/3 + 4");
+disp(x/3 + 4);
+disp("\nsin(x) + 4");
+disp(sin(x) + 4);
+disp("\nsin(x + x)");
+disp(sin(x + x));
+disp("\nsin(x + 2)");
+disp(sin(x + 2));
+disp("\nsin(2*x)");
+disp(sin(2*x));
+disp("\nsin(5*x + 3)");
+disp(sin(5*x + 3));
+disp("\nsin(cos(x))");
+disp(sin(cos(x)));
+disp("\nsin(cos(x + 2))");
+disp(cos(sin(x + 2)));
+disp("\nsin(cos(2*x + 2))");
+disp(sin(cos(2*x + 2)));
\ No newline at end of file
diff --git a/domain.m b/domain.m
new file mode 100644
index 0000000000..d6fbeaceab
--- /dev/null
+++ b/domain.m
@@ -0,0 +1,7 @@
+function out = domain(varargin)
+%DOMAIN  fake constructor for a fake DOMAIN class.
+  for i=1:nargin
+    varargin{i} = double(varargin{i});
+  end
+  out = horzcat(varargin{:});
+end
diff --git a/octave_chebfun_test_MWE.m b/octave_chebfun_test_MWE.m
new file mode 100644
index 0000000000..08bf6b5a91
--- /dev/null
+++ b/octave_chebfun_test_MWE.m
@@ -0,0 +1,42 @@
+% A container for Octave tests
+% As we try to port Chebfun to Octave, let's keep track of what works.
+% Eventually we have a goal to run the whole Chebfun test suite but small steps first.
+
+%!test
+%! % make a chebfun
+%! x = chebfun('x');
+%! assert (isa(x, 'chebfun'))
+
+%!test
+%! f = chebfun('x', [-2 2]);
+%! assert (isequal (f.domain, [-2 2]))
+
+%!test
+%! f = chebfun('x + x', [-5 5]);
+%! assert (isequal (f.pointValues, [-10;10]))
+
+%!test
+%! f = chebfun('2*x + 3', [-1 1]);
+%! assert (isequal (f(1), 5))
+%! assert (isequal (f(0), 3))
+
+%!test
+%! f = chebfun('x/3 + 4', [-10 10]);
+%! assert (isequal (f(0), 4))
+%! assert (isequal (f(3), 5))
+
+%!test
+%! f = chebfun('sin(x)', [0 pi]);
+%! area = integral(f);
+%! assert (area, 2, -2*eps)
+%assert (abs(area - 2) <= 4*eps)
+
+%!test
+%! f = chebfun('sin(x) + 2', [0 pi]);
+%! area = integral(f);
+%! assert (area, 8.283185307179588, -2*eps)
+
+%!test
+%! f = chebfun('sin(cos(2*x + 2))', [0 pi]);
+%! area = integral(f);
+%! assert (area, 0, -2*eps)
\ No newline at end of file
diff --git a/octave_notes.txt b/octave_notes.txt
new file mode 100644
index 0000000000..52d9214d9d
--- /dev/null
+++ b/octave_notes.txt
@@ -0,0 +1,42 @@
+Relevant Octave bugs
+--------------------
+
+https://savannah.gnu.org/bugs/?44035
+
+https://savannah.gnu.org/bugs/?51285
+
+Probably others!
+
+
+- - - - - - - - - - - - - - - -
+
+Older notes
+-----------
+
+Issue: class-related functions not same as private static methods
+
+http://www.mathworks.co.uk/help/matlab/matlab_oop/specifying-methods-and-functions.html?refresh=true#br2la89
+
+private static methods should use dot notation
+
+My interpretation of the Matlab docs if that calling a static method
+should require the class name (or instance) and dot notation
+("myclass.mystaticmethod()").
+
+Citations:
+
+[http://www.mathworks.co.uk/help/matlab/matlab_oop/static-methods.html]
+
+> % Calling a static method requires the class name
+
+[http://www.mathworks.co.uk/help/matlab/matlab_oop/developing-classes--typical-workflow.html]
+
+Matlab seems to allow you to be sloppy in this way inside your own class
+(or at least I can't find it documented anywhere that you are allowed to do this.)
+
+
+If really you want to use functions without the dot, there is another
+feature: "class-related function" which are also private to the m-file
+(and not e.g., methods of the class defined in other m-files).
+
+