diff --git a/README.md b/README.md
index ff076e6..11b3699 100644
--- a/README.md
+++ b/README.md
@@ -1,3 +1,298 @@
-# EllipsesAndSpheroids
+# On intersections of ellipses or spheroids with a shared focal point
 
-Cool stuff on intersections of ellipses and spheroids with a shared focal point will appear here in Feb 2024.
+Finding intersections of [ellipses](https://en.wikipedia.org/wiki/Ellipse) or [spheroids](https://en.wikipedia.org/wiki/Spheroid) is a fairly complicated task which typically involves solving a quartic equation with an iterative solver.
+However, the story is very different if the ellipses/spheroids share a focal point.
+We show below that the problem of finding the intersections reduces to the solution of a 4-by-4 (for spheroids) or a 3-by-3 (for ellipses) linear system. No iterative solver is needed. Also Matlab / Octave code for the solver and plenty of examples are given.
+
+Let us begin with a couple of pictures of ellipses with a shared focal point. The intersections are shown as black stars and the shared focal point as a black circle.
+
+![example-two-ellipses](pictures/example-ellipse-2.png)
+![example-three-ellipses](pictures/example-ellipse-3.png)
+
+## Theory
+
+We describe the theory below only for spheroids. Ellipses are just 2D special cases of spheroids. Exactly same equations hold for ellipses by dropping the Z-coodinates.
+
+Without loss of generality, we can assume that the shared focal point lies at the origin.
+Assume we have $n$ spheroids. Denote the other focal point of spheroid $k=1...n$ by $f^k = (f_x^k, f_y^k, f_y^z)$.
+Assume that $p = (x, y, z)$ is a point where all $n$ spheroids intersect.
+
+Let us denote the distance from $p$ to the origin (i.e. to the shared focal point) by $w$. That is, $w = \vert\vert p \vert\vert$.
+Furthermore, denote the distance from $p$ to the other focal point $f^k$ by $v_k$. That is, $v_{k} = \vert\vert p - f^k \vert\vert$.
+
+By definition of the spheroid, the diameter of the $k$'th spheroid is $d_k = w + v_k$, meaning that $v_k = d_k - w$.
+The leads us to the following system of $n+1$ equations.
+
+$$\qquad\begin{cases}
+ x^2 + y^2 + z^2 &= w^2 \\
+{\left ( x - f^1_x \right )}^2 + {\left ( y - f^1_y \right )}^2 + {\left ( x - f^1_z \right )}^2 &= v_1^2 = {\left ( d_1 - w \right )}^2 \\
+  & \vdots \\
+{\left ( x - f^n_x \right )}^2 + {\left ( y - f^n_y \right )}^2 + {\left ( x - f^n_z \right )}^2 &= v_n^2 = {\left ( d_n - w \right )}^2  \\
+\end{cases}\qquad\qquad\qquad (1)$$
+
+Now expand the squares and subtract the first equation from the others. This gives a system of $n$ equations as follows.
+
+$$\qquad\begin{cases}
+-2 x f^1_x  + \left ({f^1_x}\right )^2 -2 y f^1_y + \left ({f^1_y}\right )^2 -2 z f^1_z + \left ({f^1_z}\right )^2 &=  d_1^2 - 2 d_1 w  \\
+  & \vdots \\
+-2 x f^n_x  + \left ({f^n_x}\right )^2 -2 y f^n_y + \left ({f^n_y}\right )^2 -2 z f^n_z + \left ({f^n_z}\right )^2 &=  d_n^2 - 2 d_n w  \\
+\end{cases}\qquad\qquad (2)$$
+
+But this is a system of $n$ linear equations with the unknown vector being ${\left ( x \ y \ z \ w \right )}^T$ like so
+
+$$\qquad
+\left(
+ \begin{matrix}
+  -2 f^1_x & -2 f^1_y & -2 f^1_z & -2 d_1 \\
+  \vdots  & \vdots  & \vdots & \vdots  \\
+  -2 f^n_x & -2 f^n_y & -2 f^n_z & -2 d_n
+ \end{matrix}\right)
+\left(
+\begin{matrix}
+  x \\
+  y \\
+  z \\
+  w
+\end{matrix}
+\right) =
+\left(
+  \begin{matrix}
+   d_1^2 - \left ({f^1_x}\right )^2 - \left ({f^1_y}\right )^2- \left ({f^1_z}\right )^2 \\
+   \vdots \\
+   d_n^2 - \left ({f^n_x}\right )^2 - \left ({f^n_y}\right )^2- \left ({f^n_z}\right )^2 \\
+  \end{matrix}
+\right)
+\qquad\qquad (3)$$
+
+If $n < 3$ (for spheroids), there are infinitely many solutions as the intersection is a curve in 3D space.
+On the other hand 4 spheroids (or 3 ellipses) are enough to determine a unique intersection point.
+So the number of spheroids of interest is $n=3$ or 4. Likewise, the interesting number of ellipses is $n=2$ or 3.
+We will show below that for both ellipses and spheroids with a shared focal point, the number of intersections is always 0, 1 or 2.
+It follows from the fact that all intersections are either the solution of the linear problem (3) or, as we will see, they can be solved from roots of a 2nd order polynomial.
+
+If there are only 3 spheroids, meaning that the matrix on the left in equation (3) has only 3 rows, we can append a row of 4 zeros at the bottom to make it a 4-by-4 matrix.
+Likewise, we append a zero to the end of the vector on the right.
+
+Let's denote the matrix on the left-hand side of (3) by $\mathbf{T}$ and the vector on the right by $r$.
+So equation (3) can be written as $\mathbf{T}\ u = r$ where $u^T = \left ( x \ y \ z \ w \right )$.
+
+### Full rank matrix
+
+If matrix $\mathbf{T}$ has full rank, the problem is solved with an ordinary matrix inversion. That is, $u = \mathbf{T} \setminus r$ in Matlab lingo.
+Now $u(1:3) = \left [ x, y, z \right ]$ will be the intersection point and $u(4) = w$ is the distance between the intersection point and
+the shared focal point (which was assumed to be the origin). That is, $w^2 = x^2 + y^2 + z^2$.
+
+The matrix can have a full rank only if there are 4 spheroids intersect at one point only. Otherwise, the matrix will be rank deficient. We examine the solutions of the rank deficient cases below.
+
+### Rank deficient matrix
+
+Let's split matrix $\mathbf{T}$ into four 4-by-1 column vectors and denote them by $A$, $B$, $C$ and $D$. So equation (3) becomes
+
+$$\qquad\qquad\left [ A \ B \ C \ D \right ] u = r \qquad\qquad\qquad\qquad (4)$$
+
+If $\mathbf{T}$ is rank deficient and the rank is 3 (for spheroids) or 2 (for ellipses), the problem of finding the coordinates of the intersection points will reduce to finding roots of a 2nd order polynomial as we will see later. If the rank is less than that, the intersections are curves which are out of scope of this document.
+
+Let's continue with spheroids. If the rank is 3, there are 4 ways the column vectors $A, B, C, D$ can depend on each other.
+We'll go over each of these cases. But first we calculate QR-decomposition of matrix $T$ with $[\mathbf{Q}, \mathbf{R}] = qr(\mathbf{T})$.
+The rank of $\mathbf{T}$  equals the number of non-zero elements on the diagonal of the upper triangular matrix $\mathbf{R}$.
+Hence, if the rank is 3, there must be one and only one zero on the diagonal of $\mathbf{R}$.
+
+#### Case 1: $A$ is zero
+
+Assume that the zero is at the first element of the diagonal, so $\mathbf{R}_ {1,1}=0$. The means that the first column $A = {\left ( 0 \ 0 \ 0 \ 0 \right )}^T$.
+
+In this case we can forget column $A$ and the $x$-coordinate for while and solve $y$, $z$ and $w$ from the 4-by-3 linear system
+
+$$\qquad\qquad\left [ B \ C \ D \right ] {\left (y \ z \ w \right )}^T = r \qquad\qquad\qquad (5)$$
+
+This can be done with Matlab solver like so: $yzw = \mathbf{T}(:, 2:4) \setminus r$.
+
+But the first equation in (1) implies that $x = \pm \sqrt{w^2-y^2-z^2}$. This leaves us with two possible intersection points.
+
+#### Case 2: $B$ is a multiple of $A$
+
+Assume that the zero is at the second element of the diagonal, so $\mathbf{R}_ {2,2}=0$.
+This means that $B = g_A \ A$ for some real number $g_A$.
+
+As it happens, we can solve it simply as $g_A=\mathbf{R}_ {1,2} / \mathbf{R}_ {2,2}$.
+Using $B = g_A A$ in (4) we get
+
+$$\qquad\begin{matrix}
+A x + g_A A y + C z + D w &= r \\
+\Rightarrow A (x + g_A\ y) + C z + D w &= r
+\end{matrix}$$
+
+Now denote $t_A := x + g_A\ y$, meaning that $x = t_A - g_A\ y$. Then solve the 4-by-3 linear system $[A \ C \ D] \ {\left (t_A \ z \ w \right )}^T = r$ for $t_A$, $z$ and $w$.
+
+Recall from (1) that
+
+$$\qquad\begin{matrix}
+w^2 = x^2 + y^2 + z^2 = {\left ( t_A - g_A\ y \right )}^2 + y^2 + z^2 \\
+\Rightarrow  \left ( g_A^2 + 1 \right )\ y^2 - 2 g_A\ t_A\ y + \left ( t_A^2 + z^2 - w^2 \right ) = 0 \\
+\end{matrix}\qquad\qquad\qquad (6)$$
+
+This is a 2nd order polynomial of $y$ so we get two possible solutions for $y$. $z$ and $w$ were solved above and $x = t_A - g_A\ y$, which gives an $x$ for each solution of $y$.
+
+#### Case 3: $C$ is a linear combination of $A$ and $B$
+
+Assume that the zero is at the third element of the diagonal, so $\mathbf{R}_ {3,3}=0$.
+This means that $C = g_A \ A + g_B \ B$ for some real numbers $g_A$ and $g_B$.
+
+We can solve $g_A$ and $g_B$ from $\mathbf{R}$ by noting that (using Matlab notation)
+$\mathbf{R} (1:2,3)\ {\left ( g_A \ g_B \right )}^T = \mathbf{R} (1:2,1:2)$. So, we solve a 2-by-2 upper triangular system
+
+$$\qquad\qquad\left (
+\begin{matrix}
+g_A \\
+g_B
+\end{matrix}\right ) = \mathbf{R} (1:2,1:2) \setminus \mathbf{R} (1:2,3)$$
+
+Using $C = g_A \ A + g_B \ B$ in (4) we get
+
+$$\qquad\begin{matrix}
+A x + B y + C z + D w &= r \\
+\Rightarrow A (x + g_A\ z) + B (y + g_B\ z) + D w &= r
+\end{matrix}$$
+
+Now denote
+
+$$\qquad\begin{matrix}
+t_A := x + g_A\ z & \Rightarrow x = t_A - g_A\ z \\
+t_B := y + g_B\ z & \Rightarrow y = t_B - g_B\ z \\
+\end{matrix}$$
+
+Then solve the 4-by-3 linear system $[A \ B \ D] \ {\left (t_A \ t_B \ w \right )}^T = r$ for $t_A$, $t_B$ and $w$.
+
+Recall from (1) that
+
+$$\qquad\begin{matrix}
+w^2 = x^2 + y^2 + z^2 = {\left ( t_A - g_A\ z \right )}^2 + {\left ( t_B - g_B\ z \right )}^2 + z^2 \\
+\Rightarrow  \left ( g_A^2 + g_B^2 + 1 \right )\ z^2 - 2 \left (g_A\ t_A\ + g_B\ t_B\ \right )\ z + \left ( t_A^2 + t_B^2 - w^2 \right ) = 0 \\
+\end{matrix}\qquad\qquad\qquad (7)$$
+
+This is a 2nd order polynomial of $z$ so we get two possible solutions for $z$. $w$ was solved above and $x = t_A - g_A\ z$ and $y = t_B - g_B\ z$ which give a $x$ and $y$ for each solution of $z$.
+
+#### Case 4: $D$ is a linear combination of $A$, $B$ and $C$
+
+Assume that the zero is at the fourth element of the diagonal, so $\mathbf{R}_ {4,4}=0$.
+Notice that this happens typically when we have only 3 spheroids, meaning that we have added a fourth row of zeros to matrix $\mathbf{T}$ above to make it 4-by-4.
+
+This means that $D = g_A \ A + g_B \ B + g_C \ C$ for some real numbers $g_A$, $g_B$ and $g_C$.
+
+
+We can solve $g_A$, $g_B$ and $g_C$ from $\mathbf{R}$ by noting that (using Matlab notation)
+$\mathbf{R} (1:3,4)\ {\left ( g_A \ g_B \ g_C \right )}^T = \mathbf{R} (1:3,1:3)$. So, we solve a 3-by-3 upper triangular system
+
+$$\qquad\qquad\left (
+\begin{matrix}
+g_A \\
+g_B \\
+g_C
+\end{matrix}\right ) = \mathbf{R} (1:3,1:3) \setminus \mathbf{R} (1:3,4)$$
+
+Using $D = g_A \ A + g_B \ B + g_C \ C$ in (4) we get
+
+$$\qquad\begin{matrix}
+A x + B y + C z + D w &= r \\
+\Rightarrow A (x + g_A\ w) + B (y + g_B\ w) + C (z + g_C\ w) &= r
+\end{matrix}$$
+
+Now denote
+
+$$\qquad\begin{matrix}
+t_A := x + g_A\ w & \Rightarrow x = t_A - g_A\ w \\
+t_B := y + g_B\ w & \Rightarrow y = t_B - g_B\ w \\
+t_C := z + g_C\ w & \Rightarrow z = t_C - g_C\ w \\
+\end{matrix}$$
+
+Then solve the 4-by-3 linear system $[A \ B \ C] \ {\left (t_A \ t_B \ t_C \right )}^T = r$ for $t_A$, $t_B$ and $t_C$.
+
+Recall from (1) that
+
+$$\qquad\begin{matrix}
+w^2 = x^2 + y^2 + z^2 = {\left ( t_A - g_A\ w \right )}^2 + {\left ( t_B - g_B\ w \right )}^2 + {\left ( t_C - g_C\ w \right )}^2 \\
+\Rightarrow  \left ( g_A^2 + g_B^2 + g_C^2 - 1 \right )\ w^2 - 2 \left (g_A\ t_A\ + g_B\ t_B + g_C\ t_C\ \right )\ w + \left ( t_A^2 + t_B^2 + t_C^2  \right ) = 0 \\
+\end{matrix}\qquad\qquad\qquad (8)$$
+
+This is a 2nd order polynomial of $w$ so we get two possible solutions for $w$.
+Equations $x = t_A - g_A\ wy$, $y = t_B - g_B\ w$ and $z = t_C - g_C\ w$ give a $x$, $y$ and $z$ for each solution of $w$.
+
+**Note**: in every case 1...4 above, it is possible that the zeros of the 2nd order polynomials are complex numbers. If this is the case, the ellipses / spheroids do not intersect. This can happen for example if they are concentric.
+
+## Matlab / Octave solver
+
+The solver can be found in file [solveEllipseIntersections.m](octave/solveEllipseIntersections.m).
+The same solver works for both ellipses and spheroids. If the input coordinate vectors have 2 elements, the solver works in 2D mode solving intersections of ellipses. If they have 3 elements, it works in 3D mode for spheroids.
+
+The prototype of the function is
+
+````MATLAB
+    function [intersectionPoints, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters)
+````
+
+**Inputs**
+- `sharedFocalPoint` $~$  Coordinate of the focal point shared by all ellipses / spheroids. A $2\times 1$ vector for ellipses or a $3\times 1$ for spheroids.
+- `otherFocalPoints` $~$  The other $N$ (non-shared) focal points as a $2\times N$ or $3\times N$ matrix. In 2D mode, $N$ must be either 2 or 3. In 3D mode, $N$ must be either 3 or 4. So you can solve intersections of either 2 or 3 ellipses, or 3 or 4 spheroids.
+- `diameters` $~~$  Diameters of the ellipses / spheroids along the major semiaxes as an $N$-by-1 vector.
+
+**Outputs**
+- `intersectionPoints` $~$ Coordinates of the points where the $N$ ellipses / spheroids intersect as column vectors. If there are $M$ intersection points, the output is a $2\times M$ or $3\times M$ matrix. $M$ is either 0, 1 or 2. So the matrix is empty if the ellipses or spheroids do not intersect at all.
+- `relativeModelError` $~$ (_optional_) Relative difference between the actual squared distance of the calculated intersection point(s)
+                            from the shared focal point and the squared distance predicted by the model.
+                            In terms of equations, $\epsilon =  \sqrt\frac{\vert x^2+y^2+z^2 - w^2 \vert}{x^2+y^2+z^2}$ where $x,y,x$ is the relative position of the intersection point with respect to the shared focal point.
+                            If the relative error is 'large', say, bigger than 0.01, there is reason to believe that accuracy of the returned intersection points is not perfect. The reason can for instance be measurement errors in the input diameters.
+
+### First example
+
+Here is a example on how to solve intersections of two ellipses.
+
+````MATLAB
+    % Coordinates of the shared focal point as a column vector
+    sharedFocalPoint = [0 0]';
+
+    % The other focal points of two ellipses at (2,0) and (1,1) as column vectors
+    otherFocalPoints = [2 0; 1 1]';
+
+    % Diameters of the ellipses
+    diameters = [3 4]';
+
+    % Calculate intersections
+    [intersections, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters)
+
+    % Show the ellipses
+    plotEllipses(sharedFocalPoint, otherFocalPoints, diameters, intersections, "Ellipse demo");
+````
+
+The output will be
+
+````MATLAB
+intersections =
+
+   1.5817   2.4183
+  -1.0306   0.3639
+
+relativeModelError =
+
+   2.4961e-08   3.2242e-08
+````
+
+![ellipse-demo](pictures/ellipse-demo.png)
+
+### More examples
+
+There are scripts which give examples of all four rank deficient cases described in earlier chapters. Matlab and Octave scripts are living in separate folders because the way they produce graphics output is slightly different.
+
+- [demo2D_2.m](octave/demo2D_2.m) has examples on intersections of 2 ellipses. The pictures generated by this script are [2-ellipses-case-1.png](pictures/2-ellipses-case-1.png), [2-ellipses-case-2.png](pictures/2-ellipses-case-2.png) and [2-ellipses-case-3.png](pictures/2-ellipses-case-3.png).
+- [demo2D_3.m](octave/demo2D_3.m) has examples on intersections of 3 ellipses. The pictures generated by this script are [3-ellipses-case-1.png](pictures/3-ellipses-case-1.png), [3-ellipses-case-2.png](pictures/3-ellipses-case-2.png) and [3-ellipses-case-3.png](pictures/3-ellipses-case-3.png). The last picture shows the case where there is only one intersection point.
+- [demo3D_3.m](octave/demo3D_3.m) has examples on intersections of 3 spheroids. The pictures generated by this script are [3-spheroids-case-1.png](pictures/3-spheroids-case-1.png), [3-spheroids-case-2.png](pictures/3-spheroids-case-2.png), [3-spheroids-case-3.png](pictures/3-spheroids-case-3.png) and [3-spheroids-case-4.png](pictures/3-spheroids-case-4.png).
+- [demo3D_4.m](octave/demo3D_4.m) has examples on intersections of 4 spheroids. The pictures generated by this script are [4-spheroids-case-1.png](pictures/4-spheroids-case-1.png), [4-spheroids-case-2.png](pictures/4-spheroids-case-2.png), [4-spheroids-case-3.png](pictures/4-spheroids-case-3.png) and [4-spheroids-case-4.png](pictures/4-spheroids-case-4.png). The last picture shows the case where there is only one intersection point.
+
+
+## Effect of measurement errors
+
+Assume that the diameters are measured somehow and the measurements contain small errors. In this case it may happen that the ellipses / spheroids do not appear to intersect at the same point. Let's take a look at an example where noise has been added to diameters of three ellipses which otherwise intersect at a single point (like in [this picture](pictures/3-ellipses-case-3.png)).
+
+![noisy-diameters](pictures/3-ellipses-case-4.png)
+
+Despite the errors in the diameter lengths, the intersection point given by the solver is still 'near' the true intersection point. However, it is not possible to say exactly how near. All we can say is that the less error there is in the measured diameters, the closer the solution is to the true intersection point.
+We also know that the more error there is in the measurements, the larger the relative model error output is. But again, the relative model error does not say how far the given solution is from the correct one. All we can say is that the smaller the error output, the better.
diff --git a/matlab/demo2D_2.m b/matlab/demo2D_2.m
new file mode 100644
index 0000000..ee9264e
--- /dev/null
+++ b/matlab/demo2D_2.m
@@ -0,0 +1,80 @@
+%% This test demonstrates intersections of 2 ellipses
+
+% Coordinate of the shared focal point
+sharedFocalPoint = [-10; -20];
+
+% Coordinates of the other focal points of 2 ellipses.
+numEllipses = 2;
+otherFocalPoints = zeros(2, numEllipses);
+
+%% (Case 1: The major axes are collinear and upright)
+otherFocalPoints(:,1) = [0; 20] + sharedFocalPoint;
+otherFocalPoints(:,2) = [0; 10] + sharedFocalPoint;
+
+% Select point p at which the ellipses intersect and determine the diameters accordingly.
+% There will be 2 points of intersection, one of which will be p.
+p = [4; 3] + sharedFocalPoint;
+diameters = threePointDistance(otherFocalPoints, p, sharedFocalPoint)';
+
+% Solve intersections and an estimate on the model error.
+% The bigger the error, the further away the estimated intersections are from the real ones. 
+% The model error is not exacly the actual distance between 
+% the estimate and the correct intersection so it should be taken with a pinch of salt.
+[intersections, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters);
+
+assert(min(vecnorm(intersections - p)) < 1e-6, "Intersection point is not correct")
+assert(max(relativeModelError) < 1e-6, "Model error in the solver is too large")
+
+% Calculate ellipse coordinates and plot the coordinates and the intersection points.
+plotEllipses(sharedFocalPoint, otherFocalPoints, diameters, intersections, "2 ellipses, Case 1");
+
+disp("Case 1: intersection coordinates as columns: ");
+disp(intersections);
+
+%% (Case 2: The major axes are collinear and tilted)
+otherFocalPoints(:,1) = [-10; 20] + sharedFocalPoint;
+otherFocalPoints(:,2) = [-5; 10] + sharedFocalPoint;
+
+% Select point p at which the ellipses intersect and determine the diameters accordingly.
+% There will be 2 points of intersection, one of which will be p.
+p = [4; 3] + sharedFocalPoint;
+diameters = threePointDistance(otherFocalPoints, p, sharedFocalPoint)';
+
+% Solve intersections and an estimate on the model error.
+% The bigger the error, the further away the estimated intersections are from the real ones. 
+% The model error is not exacly the actual distance between 
+% the estimate and the correct intersection so it should be taken with a pinch of salt.
+[intersections, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters);
+
+assert(min(vecnorm(intersections - p)) < 1e-6, "Intersection point is not correct")
+assert(max(relativeModelError) < 1e-6, "Model error in the solver is too large")
+
+% Calculate ellipse coordinates and plot the coordinates and the intersection points.
+plotEllipses(sharedFocalPoint, otherFocalPoints, diameters, intersections, "2 ellipses, Case 2");
+
+disp("Case 2: intersection coordinates as columns:");
+disp(intersections);
+
+%% (Case 3: General case with 2 ellipses)
+otherFocalPoints(:,1) = [10; 4] + sharedFocalPoint;
+otherFocalPoints(:,2) = [0; 10] + sharedFocalPoint;
+
+% Select point p at which the ellipses intersect and determine the diameters accordingly.
+% There will be 2 points of intersection, one of which will be p.
+p = [4; 3] + sharedFocalPoint;
+diameters = threePointDistance(otherFocalPoints, p, sharedFocalPoint)';
+
+% Solve intersections and an estimate on the model error.
+% The bigger the error, the further away the estimated intersections are from the real ones. 
+% The model error is not exacly the actual distance between 
+% the estimate and the correct intersection so it should be taken with a pinch of salt.
+[intersections, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters);
+
+assert(min(vecnorm(intersections - p)) < 1e-6, "Intersection point is not correct")
+assert(max(relativeModelError) < 1e-6, "Model error in the solver is too large")
+
+% Calculate ellipse coordinates and plot the coordinates and the intersection points.
+plotEllipses(sharedFocalPoint, otherFocalPoints, diameters, intersections, "2 ellipses, Case 3");
+
+disp("Case 3: intersection coordinates as columns:");
+disp(intersections);
\ No newline at end of file
diff --git a/matlab/demo2D_3.m b/matlab/demo2D_3.m
new file mode 100644
index 0000000..9857215
--- /dev/null
+++ b/matlab/demo2D_3.m
@@ -0,0 +1,112 @@
+%% This test demonstrates intersections of 3 ellipses
+
+% Coordinate of the shared focal point
+sharedFocalPoint = [-10; -20];
+
+% Coordinates of the other focal points of 2 ellipses.
+numEllipses = 3;
+otherFocalPoints = zeros(2, numEllipses);
+
+%% (Case 1: The major axes are collinear and upright)
+otherFocalPoints(:,1) = [0; 20] + sharedFocalPoint;
+otherFocalPoints(:,2) = [0; 10] + sharedFocalPoint;
+otherFocalPoints(:,3) = [0; -16] + sharedFocalPoint;
+
+% Select point p at which the ellipses intersect and determine the diameters accordingly.
+% There will be 2 points of intersection, one of which will be p.
+p = [4; 3] + sharedFocalPoint;
+diameters = threePointDistance(otherFocalPoints, p, sharedFocalPoint)';
+
+% Solve intersections and an estimate on the model error.
+% The bigger the error, the further away the estimated intersections are from the real ones.
+% The model error is not exacly the actual distance between
+% the estimate and the correct intersection so it should be taken with a pinch of salt.
+[intersections, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters);
+
+assert(min(vecnorm(intersections - p)) < 1e-6, "Intersection point is not correct")
+assert(max(relativeModelError) < 1e-6, "Model error in the solver is too large")
+
+% Calculate ellipse coordinates and plot the coordinates and the intersection points.
+plotEllipses(sharedFocalPoint, otherFocalPoints, diameters, intersections, "3 ellipses, Case 1");
+
+disp("Case 1: intersection coordinates as columns: ");
+disp(intersections);
+
+%% (Case 2: The major axes are collinear and tilted)
+otherFocalPoints(:,1) = [-10; 20] + sharedFocalPoint;
+otherFocalPoints(:,2) = [-5; 10] + sharedFocalPoint;
+otherFocalPoints(:,3) = [-8; 16] + sharedFocalPoint;
+
+% Select point p at which the ellipses intersect and determine the diameters accordingly.
+% There will be 2 points of intersection, one of which will be p.
+p = [4; 3] + sharedFocalPoint;
+diameters = threePointDistance(otherFocalPoints, p, sharedFocalPoint)';
+
+% Solve intersections and an estimate on the model error.
+% The bigger the error, the further away the estimated intersections are from the real ones.
+% The model error is not exacly the actual distance between
+% the estimate and the correct intersection so it should be taken with a pinch of salt.
+[intersections, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters);
+
+assert(min(vecnorm(intersections - p)) < 1e-6, "Intersection point is not correct")
+assert(max(relativeModelError) < 1e-6, "Model error in the solver is too large")
+
+% Calculate ellipse coordinates and plot the coordinates and the intersection points.
+plotEllipses(sharedFocalPoint, otherFocalPoints, diameters, intersections, "3 ellipses, Case 2");
+
+disp("Case 2: intersection coordinates as columns:");
+disp(intersections);
+
+%% (Case 3: General case with 3 ellipses)
+otherFocalPoints(:,1) = [10; 4] + sharedFocalPoint;
+otherFocalPoints(:,2) = [0; 10] + sharedFocalPoint;
+otherFocalPoints(:,3) = [10; 0] + sharedFocalPoint;
+
+% Select point p at which the ellipses intersect and determine the diameters accordingly.
+% There will be 2 points of intersection, one of which will be p.
+p = [4; 3] + sharedFocalPoint;
+diameters = threePointDistance(otherFocalPoints, p, sharedFocalPoint)';
+
+% Solve intersections and an estimate on the model error.
+% The bigger the error, the further away the estimated intersections are from the real ones.
+% The model error is not exacly the actual distance between
+% the estimate and the correct intersection so it should be taken with a pinch of salt.
+[intersections, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters);
+
+assert(min(vecnorm(intersections - p)) < 1e-6, "Intersection point is not correct")
+assert(max(relativeModelError) < 1e-6, "Model error in the solver is too large")
+
+% Calculate ellipse coordinates and plot the coordinates and the intersection points.
+plotEllipses(sharedFocalPoint, otherFocalPoints, diameters, intersections, "3 ellipses, Case 3");
+
+disp("Case 3: (One intersection only) Intersection coordinates as columns:");
+disp(intersections);
+
+%% (Case 4: Non-solvable case with 3 ellipses)
+% This is an example of the solver giving a wrong answer and returning a high
+% model error when the requirement that the ellipses intersect at one point does not hold.
+otherFocalPoints(:,1) = [10; 4] + sharedFocalPoint;
+otherFocalPoints(:,2) = [0; 10] + sharedFocalPoint;
+otherFocalPoints(:,3) = [10; 0] + sharedFocalPoint;
+
+% Select point p at which the ellipses intersect and determine the diameters accordingly.
+% There will be 2 points of intersection, one of which will be p.
+p = [4; 3] + sharedFocalPoint;
+diameters = threePointDistance(otherFocalPoints, p, sharedFocalPoint)';
+
+% The 3rd ellipse does not intersect at the same point where the 1st and 2nd ellipses intersect
+pp = p + [1; 1];
+diameters(1) = threePointDistance(otherFocalPoints(:,1), pp, sharedFocalPoint)';
+
+% Solve intersections and an estimate on the model error.
+% The bigger the error, the further away the estimated intersections are from the real ones.
+% The model error is not exacly the actual distance between
+% the estimate and the correct intersection so it should be taken with a pinch of salt.
+[intersections, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters);
+
+% Calculate ellipse coordinates and plot the coordinates and the intersection points.
+plotEllipses(sharedFocalPoint, otherFocalPoints, diameters, intersections, ...
+             ["3 ellipses intersecting at different points, rel.model error=" num2str(relativeModelError)]);
+
+disp(["Case 4: intersection with model error " num2str(relativeModelError)]);
+disp(intersections);
diff --git a/matlab/demo3D_3.m b/matlab/demo3D_3.m
new file mode 100644
index 0000000..d2525f1
--- /dev/null
+++ b/matlab/demo3D_3.m
@@ -0,0 +1,107 @@
+%% This test demonstrates intersections of 3 spheroids
+
+% Coordinate of the shared focal point
+sharedFocalPoint = [-10; -20; 30];
+
+% Coordinates of the other focal points of 2 spheroids.
+numSpheroids = 3;
+otherFocalPoints = zeros(3, numSpheroids);
+
+%% (Case 1: The major axes are on the same yz-plane where x is constant)
+otherFocalPoints(:,1) = [0; 20; 1] + sharedFocalPoint;
+otherFocalPoints(:,2) = [0; 10; -2] + sharedFocalPoint;
+otherFocalPoints(:,3) = [0; -16; 3] + sharedFocalPoint;
+
+% Select point p at which the spheroids intersect and determine the diameters accordingly.
+% There will be 2 points of intersection, one of which will be p.
+p = [4; 3; 11] + sharedFocalPoint;
+diameters = threePointDistance(otherFocalPoints, p, sharedFocalPoint)';
+
+% Solve intersections and an estimate on the model error.
+% The bigger the error, the further away the estimated intersections are from the real ones. 
+% The model error is not exacly the actual distance between 
+% the estimate and the correct intersection so it should be taken with a pinch of salt.
+[intersections, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters);
+
+assert(min(vecnorm(intersections - p)) < 1e-6, "Intersection point is not correct")
+assert(max(relativeModelError) < 1e-6, "Model error in the solver is too large")
+
+% Calculate spheroid coordinates and plot the coordinates and the intersection points.
+plotSpheroids(sharedFocalPoint, otherFocalPoints, diameters, intersections, "3 spheroids, Case 1");
+
+disp("Case 1: intersection coordinates as columns: ");
+disp(intersections);
+
+%% (Case 2: The major axes are on the same tilted plane (x and y coordinates are linearly dependent))
+otherFocalPoints(:,1) = [10; 20; 1] + sharedFocalPoint;
+otherFocalPoints(:,2) = [5; 10; -2] + sharedFocalPoint;
+otherFocalPoints(:,3) = [-8; -16; 3] + sharedFocalPoint;
+
+% Select point p at which the spheroids intersect and determine the diameters accordingly.
+% There will be 2 points of intersection, one of which will be p.
+p = [4; 3; 11] + sharedFocalPoint;
+diameters = threePointDistance(otherFocalPoints, p, sharedFocalPoint)';
+
+% Solve intersections and an estimate on the model error.
+% The bigger the error, the further away the estimated intersections are from the real ones. 
+% The model error is not exacly the actual distance between 
+% the estimate and the correct intersection so it should be taken with a pinch of salt.
+[intersections, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters);
+
+assert(min(vecnorm(intersections - p)) < 1e-6, "Intersection point is not correct")
+assert(max(relativeModelError) < 1e-6, "Model error in the solver is too large")
+
+% Calculate spheroid coordinates and plot the coordinates and the intersection points.
+plotSpheroids(sharedFocalPoint, otherFocalPoints, diameters, intersections, "3 spheroids, Case 2");
+
+disp("Case 2: intersection coordinates as columns:");
+disp(intersections);
+
+%% (Case 3: z coordinates of the major axes are linearly dependent on x- and y-coordinates.)
+otherFocalPoints(:,1) =[10; 1; 10+1] + sharedFocalPoint;
+otherFocalPoints(:,2) = [5; -2; 5-2] + sharedFocalPoint;
+otherFocalPoints(:,3) = [-8; 3; -8+3] + sharedFocalPoint;
+
+% Select point p at which the spheroids intersect and determine the diameters accordingly.
+% There will be 2 points of intersection, one of which will be p.
+p = [4; 3; 11] + sharedFocalPoint;
+diameters = threePointDistance(otherFocalPoints, p, sharedFocalPoint)';
+
+% Solve intersections and an estimate on the model error.
+% The bigger the error, the further away the estimated intersections are from the real ones. 
+% The model error is not exacly the actual distance between 
+% the estimate and the correct intersection so it should be taken with a pinch of salt.
+[intersections, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters);
+
+assert(min(vecnorm(intersections - p)) < 1e-6, "Intersection point is not correct")
+assert(max(relativeModelError) < 1e-6, "Model error in the solver is too large")
+
+% Calculate spheroid coordinates and plot the coordinates and the intersection points.
+plotSpheroids(sharedFocalPoint, otherFocalPoints, diameters, intersections, "3 spheroids, Case 3");
+
+disp("Case 3: intersection coordinates as columns:");
+disp(intersections);
+%% (Case 4: General case with 3 spheroids)
+otherFocalPoints(:,1) = [-2; -5; 1] + sharedFocalPoint;
+otherFocalPoints(:,2) = [0; 10; 2] + sharedFocalPoint;
+otherFocalPoints(:,3) = [10; 0; 3] + sharedFocalPoint;
+
+% Select point p at which the spheroids intersect and determine the diameters accordingly.
+% There will be 2 points of intersection, one of which will be p.
+p = [4; 3; 11] + sharedFocalPoint;
+diameters = threePointDistance(otherFocalPoints, p, sharedFocalPoint)';
+
+% Solve intersections and an estimate on the model error.
+% The bigger the error, the further away the estimated intersections are from the real ones. 
+% The model error is not exacly the actual distance between 
+% the estimate and the correct intersection so it should be taken with a pinch of salt.
+[intersections, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters);
+
+assert(min(vecnorm(intersections - p)) < 1e-6, "Intersection point is not correct")
+assert(max(relativeModelError) < 1e-6, "Model error in the solver is too large")
+
+% Calculate spheroid coordinates and plot the coordinates and the intersection points.
+plotSpheroids(sharedFocalPoint, otherFocalPoints, diameters, intersections, "3 spheroids, Case 4");
+
+disp("Case 4: intersection coordinates as columns:");
+disp(intersections);
\ No newline at end of file
diff --git a/matlab/demo3D_4.m b/matlab/demo3D_4.m
new file mode 100644
index 0000000..e58f727
--- /dev/null
+++ b/matlab/demo3D_4.m
@@ -0,0 +1,111 @@
+%% This test demonstrates intersections of 4 spheroids
+
+% Coordinate of the shared focal point
+sharedFocalPoint = [-10; -20; 30];
+
+% Coordinates of the other focal points of 2 spheroids.
+numSpheroids = 4;
+otherFocalPoints = zeros(3, numSpheroids);
+
+%% (Case 1: The major axes are on the same yz-plane where x is constant)
+otherFocalPoints(:,1) = [0; 20; 1] + sharedFocalPoint;
+otherFocalPoints(:,2) = [0; 10; -2] + sharedFocalPoint;
+otherFocalPoints(:,3) = [0; -16; 3] + sharedFocalPoint;
+otherFocalPoints(:,4) = [0; -14; -4] + sharedFocalPoint;
+
+% Select point p at which the spheroids intersect and determine the diameters accordingly.
+% There will be 2 points of intersection, one of which will be p.
+p = [4; 3; 11] + sharedFocalPoint;
+diameters = threePointDistance(otherFocalPoints, p, sharedFocalPoint)';
+
+% Solve intersections and an estimate on the model error.
+% The bigger the error, the further away the estimated intersections are from the real ones. 
+% The model error is not exacly the actual distance between 
+% the estimate and the correct intersection so it should be taken with a pinch of salt.
+[intersections, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters);
+
+assert(min(vecnorm(intersections - p)) < 1e-6, "Intersection point is not correct")
+assert(max(relativeModelError) < 1e-6, "Model error in the solver is too large")
+
+% Calculate spheroid coordinates and plot the coordinates and the intersection points.
+plotSpheroids(sharedFocalPoint, otherFocalPoints, diameters, intersections, "4 spheroids, Case 1");
+
+disp("Case 1: intersection coordinates as columns: ");
+disp(intersections);
+
+%% (Case 2: The major axes are on the same tilted plane (x and y coordinates are linearly dependent))
+otherFocalPoints(:,1) = [10; 20; 1] + sharedFocalPoint;
+otherFocalPoints(:,2) = [5; 10; -2] + sharedFocalPoint;
+otherFocalPoints(:,3) = [-8; -16; 3] + sharedFocalPoint;
+otherFocalPoints(:,4) = [-7; -14; -4] + sharedFocalPoint;
+
+% Select point p at which the spheroids intersect and determine the diameters accordingly.
+% There will be 2 points of intersection, one of which will be p.
+p = [4; 3; 11] + sharedFocalPoint;
+diameters = threePointDistance(otherFocalPoints, p, sharedFocalPoint)';
+
+% Solve intersections and an estimate on the model error.
+% The bigger the error, the further away the estimated intersections are from the real ones. 
+% The model error is not exacly the actual distance between 
+% the estimate and the correct intersection so it should be taken with a pinch of salt.
+[intersections, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters);
+
+assert(min(vecnorm(intersections - p)) < 1e-6, "Intersection point is not correct")
+assert(max(relativeModelError) < 1e-6, "Model error in the solver is too large")
+
+% Calculate spheroid coordinates and plot the coordinates and the intersection points.
+plotSpheroids(sharedFocalPoint, otherFocalPoints, diameters, intersections, "4 spheroids, Case 2");
+
+disp("Case 2: intersection coordinates as columns:");
+disp(intersections);
+
+%% (Case 3: z coordinates of the major axes are linearly dependent on x- and y-coordinates.)
+otherFocalPoints(:,1) =[10; 1; 10+1] + sharedFocalPoint;
+otherFocalPoints(:,2) = [5; -2; 5-2] + sharedFocalPoint;
+otherFocalPoints(:,3) = [-8; 3; -8+3] + sharedFocalPoint;
+otherFocalPoints(:,4) = [-7; -4; -7-4] + sharedFocalPoint;
+
+% Select point p at which the spheroids intersect and determine the diameters accordingly.
+% There will be 2 points of intersection, one of which will be p.
+p = [4; 3; 11] + sharedFocalPoint;
+diameters = threePointDistance(otherFocalPoints, p, sharedFocalPoint)';
+
+% Solve intersections and an estimate on the model error.
+% The bigger the error, the further away the estimated intersections are from the real ones. 
+% The model error is not exacly the actual distance between 
+% the estimate and the correct intersection so it should be taken with a pinch of salt.
+[intersections, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters);
+
+assert(min(vecnorm(intersections - p)) < 1e-6, "Intersection point is not correct")
+assert(max(relativeModelError) < 1e-6, "Model error in the solver is too large")
+
+% Calculate spheroid coordinates and plot the coordinates and the intersection points.
+plotSpheroids(sharedFocalPoint, otherFocalPoints, diameters, intersections, "4 spheroids, Case 3");
+
+disp("Case 3: intersection coordinates as columns:");
+disp(intersections);
+%% (Case 4: General case with 3 spheroids)
+otherFocalPoints(:,1) = [-2; -5; 1] + sharedFocalPoint;
+otherFocalPoints(:,2) = [0; 10; 2] + sharedFocalPoint;
+otherFocalPoints(:,3) = [10; 0; 3] + sharedFocalPoint;
+otherFocalPoints(:,4) = [-8; -3; 4] + sharedFocalPoint;
+
+% Select point p at which the spheroids intersect and determine the diameters accordingly.
+% There will be 2 points of intersection, one of which will be p.
+p = [4; 3; 11] + sharedFocalPoint;
+diameters = threePointDistance(otherFocalPoints, p, sharedFocalPoint)';
+
+% Solve intersections and an estimate on the model error.
+% The bigger the error, the further away the estimated intersections are from the real ones. 
+% The model error is not exacly the actual distance between 
+% the estimate and the correct intersection so it should be taken with a pinch of salt.
+[intersections, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters);
+
+assert(min(vecnorm(intersections - p)) < 1e-6, "Intersection point is not correct")
+assert(max(relativeModelError) < 1e-6, "Model error in the solver is too large")
+
+% Calculate spheroid coordinates and plot the coordinates and the intersection points.
+plotSpheroids(sharedFocalPoint, otherFocalPoints, diameters, intersections, "4 spheroids, Case 4");
+
+disp("Case 4: intersection coordinates as columns:");
+disp(intersections);
\ No newline at end of file
diff --git a/matlab/plotEllipses.m b/matlab/plotEllipses.m
new file mode 100644
index 0000000..f9d9c39
--- /dev/null
+++ b/matlab/plotEllipses.m
@@ -0,0 +1,84 @@
+function output = plotEllipses(sharedFocalPoint, otherFocalPoints, diameters, intersections, titleText)
+%plotEllipses Plot ellipses with the given focal points and diameters. Also show the points of intersections.
+% Inputs:
+%   sharedFocalPoint - Coordinate of shared focal points of every ellipse (2x1)
+%   otherFocalPoints - Coordinates of N other focal points (2xN) 
+%   diameters - Diameters of the ellipses along the major axis (Nx1)
+%   intersections - Coordinates of M intersections (2xM)
+
+    if nargin < 5
+        titleText = "";
+    end
+    
+    numEllipses = size(otherFocalPoints, 2);
+    numPoints = 100;
+    ellipseCoord = zeros(2,numPoints, numEllipses);
+    theta = linspace(0, 2*pi, numPoints);
+    diffCoord = otherFocalPoints - sharedFocalPoint;
+    for ii = 1:numEllipses
+        a = diameters(ii)/2;  % Major axis
+        b = sqrt(a*a - norm(diffCoord(:,ii)/2)^2); % Minor axis
+        phase = diffCoord(1,ii) + 1i*diffCoord(2,ii);
+        phase = phase / norm(phase); % Rotation of the ellipse
+        center = (otherFocalPoints(:,ii) + sharedFocalPoint) / 2; % Center of the ellipse
+        for jj = 1:numPoints
+            coord = (a * cos(theta(jj))) + 1i*(b * sin(theta(jj))); % axis aligned
+            coord = coord * phase;  % rotate
+            coord = coord + (center(1) + 1i*center(2));
+            ellipseCoord(:,jj,ii) = [real(coord); imag(coord)];
+        end
+    end
+    rgb=['g','r','b'];
+    limX = [floor(min(ellipseCoord(1,:))) ceil(max(ellipseCoord(1,:)))];
+    limY = [floor(min(ellipseCoord(2,:))) ceil(max(ellipseCoord(2,:)))];
+    legendArray = {};
+    plots = [];
+
+    figure;
+    hold on
+    % Plot the shared focal point
+    plots(end+1) = plot(sharedFocalPoint(1),sharedFocalPoint(2),'Marker','o','MarkerSize',10,'Color','k');
+    legendArray{end+1} = "Shared focal point";
+
+    % Plot non-shared focal points
+    for ii=1:numEllipses
+        plots(end+1) = plot(otherFocalPoints(1,ii),otherFocalPoints(2,ii),'Marker','*','MarkerSize',10,'Color',rgb(mod(ii,3)+1));
+        legendArray{end+1} = "focal point #" + num2str(ii);
+    end
+
+    % Plot ellipses
+    for ii=1:numEllipses
+        plots(end+1) = plot(ellipseCoord(1,:,ii), ellipseCoord(2,:,ii),'-','Color',rgb(mod(ii,3)+1));
+        legendArray{end+1} = "Ellipse #" + num2str(ii);
+    end
+
+
+    % Plot the intersections
+    for jj=1:size(intersections,2)
+        if mod(jj,2) == 0
+            mkr = 'pentagram';
+        else
+            mkr = 'hexagram';
+        end
+        plots(end+1) = plot(intersections(1,jj), intersections(2,jj),'Marker',mkr,'MarkerSize',12,'Color','k','LineStyle',':');
+        legendArray{end+1} = "Intersection #" + num2str(jj);
+    end
+
+    for jj=1:size(intersections,2)
+        for ii=1:numEllipses
+            plot([intersections(1,jj) otherFocalPoints(1,ii)], [intersections(2,jj) otherFocalPoints(2,ii)],':','Color',rgb(mod(ii,3)+1))
+            plot([intersections(1,jj) sharedFocalPoint(1)], [intersections(2,jj) sharedFocalPoint(2)],':','Color',rgb(mod(ii,3)+1))
+        end
+    end 
+
+    grid on
+    axis equal
+    ax = gca;
+    ax.XTick = limX(1):limX(2);
+    ax.YTick = limY(1):limY(2);
+    title(titleText)
+    legend(plots,legendArray);
+
+    output = gca;
+    
+end
\ No newline at end of file
diff --git a/matlab/plotSpheroids.m b/matlab/plotSpheroids.m
new file mode 100644
index 0000000..ed98164
--- /dev/null
+++ b/matlab/plotSpheroids.m
@@ -0,0 +1,97 @@
+function output = plotSpheroids(sharedFocalPoint, otherFocalPoints, diameters, intersections, titleText)
+%plotEllipses Plot spheroids with the given focal points and diameters. Also show the points of intersections.
+% Inputs:
+%   sharedFocalPoint - Coordinate of shared focal points of every ellipse (2x1)
+%   otherFocalPoints - Coordinates of N other focal points (2xN) 
+%   diameters - Diameters of the ellipses along the major axis (Nx1)
+%   intersections - Coordinates of M intersections (2xM)
+
+    if nargin < 5
+        titleText = "";
+    end
+  
+    numPoints = 50;
+    numSpheroids = size(otherFocalPoints, 2);
+    spheroidCoord = zeros(3,numPoints, numSpheroids);
+    theta = linspace(0, 2*pi, numPoints);
+    diffCoord = otherFocalPoints - sharedFocalPoint;
+
+    figure;
+    hold on
+    grid on
+    axis equal
+    legendArray = {};
+    plots = [];
+
+    for ii = 1:numSpheroids
+        a = diameters(ii)/2;  % Major axis
+        b = sqrt(a*a - norm(diffCoord(:,ii)/2)^2); % Minor axis
+        helperEllipse = zeros(3, numPoints*numPoints);
+        helperIndex = 0;
+        for jj = 1:numPoints
+            coord = [a * cos(theta(jj)) , b * sin(theta(jj))];
+            for kk = 1:numPoints
+                helperIndex = helperIndex + 1;
+                helperEllipse(:,helperIndex) = [coord(1); cos(theta(kk))*coord(2); sin(theta(kk))*coord(2)];
+            end
+        end
+
+        center = (otherFocalPoints(:,ii) + sharedFocalPoint) / 2; % Center of the spheroid
+        
+        % Cartesian x axis will be mapped into this vector,
+        % which is the direction of the major axis.
+        Vaxis = diffCoord(:,ii);
+        Vaxis = Vaxis / norm(Vaxis);
+        % Decide another axis orthogonal the Xaxis
+        [~,minIndex]=min(abs(axis));
+        switch minIndex
+            case 1
+                Uaxis = [1;0;0];
+            case 2
+                Uaxis = [0;1;0];
+            case 3
+                Uaxis = [0;0;1];
+        end
+        Waxis = Uaxis - (Uaxis' * Vaxis) * Vaxis;
+        Waxis = Waxis / norm(Waxis);
+        Uaxis=cross(Waxis, Vaxis);
+        % This matrix rotates the helper ellipse into place
+        Rot = [Vaxis Waxis Uaxis];
+    
+        spheroid = Rot * helperEllipse + center;
+    
+        plots(end+1) = plot3(spheroid(1,:), spheroid(2,:), spheroid(3,:),'.','MarkerSize',10);
+        legendArray{end+1} = "Spheroid "+num2str(ii);
+        % Plot Rx->Tx pairs
+        plots(end+1)=plot3([sharedFocalPoint(1), otherFocalPoints(1,ii)], ...
+            [sharedFocalPoint(2), otherFocalPoints(2,ii)], ...
+            [sharedFocalPoint(3), otherFocalPoints(3,ii)],'x-','LineWidth',3);
+        legendArray{end+1} = "Foci axis "+num2str(ii);
+    end
+
+    % Plot the intersections
+    numIntersections = size(intersections,2);
+    if numIntersections == 1
+            plots(end+1) = plot3(intersections(1), intersections(2), intersections(3), ...
+                'ro', 'MarkerSize',16,'MarkerFaceColor','auto','LineWidth',2);
+            legendArray{end+1} = "Intersection";
+    elseif numIntersections > 1
+        for jj=1:size(intersections,2)
+            if mod(jj,2) == 0
+                mkr = 'pentagram';
+            else
+                mkr = 'hexagram';
+            end
+            plots(end+1) = plot3(intersections(1,jj), intersections(2,jj), intersections(3,jj), ...
+                'Color','r','Marker',mkr, 'MarkerSize',16,'MarkerFaceColor','auto','LineWidth',2);
+            legendArray{end+1} = "Intersection #" + num2str(jj);
+        end
+    end
+
+    xlabel("XX")
+    ylabel("YY")
+    zlabel("ZZ")
+    legend(legendArray)
+    title(titleText)
+    output = gca; 
+end
\ No newline at end of file
diff --git a/matlab/solveEllipseIntersections.m b/matlab/solveEllipseIntersections.m
new file mode 100644
index 0000000..4c955f5
--- /dev/null
+++ b/matlab/solveEllipseIntersections.m
@@ -0,0 +1,335 @@
+function [intersectionPoints, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters)
+%solveIntersections Calculates intersections of ellipses and spheroids with a shared focal point.
+%   Calculates intersections of ellipses (2D) or spheroids (3D). All ellipses are assumed to have the first focal point  
+%   at the same position and the second focal point at the given coordinates.
+%   The algorithm for 3D spheroids is described below. The 2D case is similar but simpler.
+%
+%   The diameters along the major semiaxes determine the sum of distances from the first focal point to 
+%   any point on the locus of the ellipse or spheroid and further on to the second focal point.
+%   Without loss of generality we assume that the first focal point of each spheroid is at the origin 
+%   and fp = otherFocalPoints - sharedFocalPoint is the set of the second focal points.
+%   Now assume that [x y z] is an intersection of the spheroids. Its distance from the origin is denoted by w.
+%   Now the intersection point must satisfy the following N+1 conditions where N is the number of intersecting spheroids:
+%
+%                                           x^2 + y^2 + z^2 = w^2
+%       (x - fp(1,1))^2 + (y - fp(2,1))^2 + (z - fp(3,1))^2 = (diam(1) - w)^2
+%       ...
+%       (x - fp(1,N))^2 + (y - fp(2,N))^2 + (z - fp(3,N))^2 = (diam(N) - w)^2
+%
+% By expanding the squares in N last equations and subtracting the first equation
+% from the expanded squares, we get an N-by-4 linear system  T*[x y z w]' = r as follows:
+%
+%   [-2 fp(1,1) -2 fp(2,1) -2 fp(3,1) diam(1)] [x] = [diam(1)^2 - fp(1,1)^2 - fp(2,1)^2 - fp(3,1)^2]
+%   [                   ...                  ] [y] = [      ....                       ]
+%   [                   ...                  ] [z] = [      ....                       ]
+%   [-2 fp(1,N) -2 fp(2,N) -2 fp(3,N) diam(N)] [w] = [diam(N)^2 - fp(1,N)^2 - fp(2,N)^2 - fp(3,N)^2]
+%
+%   So T = -2 * [fp(1,:)' fp(2,:)' fp(3,:)' diam] and r =  diam.^2 - fp(1,:)'.^2 - fp(2,:)'.^2 - fp(3,:)'.^2;
+%
+%   Notice that a row of T may be linearly dependent on the other rows. Those cases can be dealt with
+%   by writing the dependent row as a linear combination of the other rows and using the fact that
+%   x^2 + y^2 + z^2 = w^2 to solve the coefficients of the linear combination. This will be explained in
+%   comments in the code below.
+%
+%   Inputs:
+%       sharedFocalPoint - The focal point which is shared by all ellipses / spheroids, DIM-by-1 vector
+%       otherFocalPoints - The other (non-shared) focal points as DIM-by-N matrix.
+%          In 2D case, N should be 2 or 3. In 3D case, N should be 3 or 4.
+%       diameters - Diameters of the ellipses along the major semiaxes as an N-by-1 vector.
+%   Output:
+%       intersectionPoints - Coordinates of the points where all N ellipses intersect. 
+%                            A 2-by-M or 3-by-M matrix where M is 1 of the intersection point is unique and 2 if there are two candidates.
+%                            If the matrix is empty, there is no common intersection point.
+%                            Note that if there is noise or other inaccuracy in the set of diameters, the all ellipses/spheroids
+%                            may not actually intersect at the given points. In this case distanceError will be large
+%                            so the inaccuracy can be detected.
+%       relativeModelError - Relative difference between the actual squared distance of the calculated intersection point(s)
+%                            from the shared focal point (i.e. x^2+y^2+z^2) and the squared distance predicted by
+%                            the model (i.e. w^2). So relativeModelError = sqrt(abs(x^2+y^2+z^2 - w^2) / x^2+y^2+z^2).
+%                            The larger the value, the less reliable the result is. 
+%                            The value is very near zero if all ellipses/spheroids intersect at the same points.
+
+assert(nargin == 3)
+% 2D or 3D mode?
+
+dimension = size(otherFocalPoints, 1);
+numFoci = size(otherFocalPoints, 2);
+
+assert(numFoci == length(diameters), "The number of focal points and the number of diameters must be the same.");
+assert(size(diameters, 2) == 1, "Diameters must be an N-by-1 vector")
+assert(size(sharedFocalPoint,1) == dimension)
+assert(size(sharedFocalPoint,2) == 1)
+
+% Normalize the coordinates so that the shared focal point is at the origin.
+focalPoints = otherFocalPoints - sharedFocalPoint;
+switch dimension
+    case 2
+        assert(numFoci == 2 || numFoci == 3, "In 2D case the number of points should be 2 or 3")
+         intersectionPoints = solveEllipseIntersections2D(focalPoints, diameters);
+    case 3
+        assert(numFoci == 3 || numFoci == 4, "In 3D case the number of points should be 3 or 4")
+        intersectionPoints = solveSpheroidIntersections3D(focalPoints, diameters);
+    otherwise
+        error("The dimension of focal points must be 2 or 3")
+end
+
+
+if isempty(intersectionPoints)
+    relativeModelError = Inf;
+    return
+end
+
+% Sometimes the imaginary components can be non-zero due to numerical errors in the solver.
+if ~isreal(intersectionPoints)
+    % The points are conjugates, so only one real part is enough.
+    intersectionPoints = real(intersectionPoints(:,1));
+end
+if nargout > 1
+    w = intersectionPoints(dimension+1,:);
+    xyz = intersectionPoints(1:dimension,:);
+    relativeModelError =  sqrt(abs(vecnorm(xyz).^2 - w.^2)) ./ (vecnorm(xyz) + realmin);
+end
+
+% Undo the normalization
+intersectionPoints = intersectionPoints(1:dimension,:) + sharedFocalPoint;
+end
+
+
+function intersectionPoints = solveEllipseIntersections2D(fp, diam)
+%   Inputs:
+%       fp - Second focal points as a 2-by-N matrix (the first focal points are always at the origin) 
+%       diam - Diameters of the ellipses along the major semiaxes, N-by-1 vector.
+%   Output:
+%       intersectionPoints - Coordinates of the points where all N ellipses intersect, or an empty vector if they don't intersect.
+%
+    DIM = 2;
+    intersectionPoints = [];
+    numFoci = length(diam); % Number of ellipses (must be DIM or DIM+1)
+
+%  An intersection point (x,y) of the spheroids and it's distance w = sqrt(x^2 + y^2) from the origin 
+% (i.e. the second focal point) satifies the following conditions:
+%       (x - fp(1,1))^2 + (y - fp(2,1))^2 = (diam(1) - w)^2
+%       ...
+%       (x - fp(1,N))^2 + (y - fp(2,N))^2 = (diam(N) - w)^2
+%                               x^2 + y^2 = w^2
+% By expanding the squares in N first equations and using the last equation, 
+% we get a linear problem T * [x y w]' = r  with matrix T and vector r defined as as follows:
+    T = -2 * [fp(1,:)' fp(2,:)' diam];
+    r = diam.^2 - fp(1,:)'.^2 - fp(2,:)'.^2;
+    
+    % QR transform will later help with sorting out corner cases later
+    [~,R] = qr(T);
+    
+    % Are there tiny elements on the diagonal of R? 
+    % If so, it means that some columns of A are linearly dependent.
+    zeroDiagIndex = findSmallIndices(diag(R));
+
+    % If there are fewer than DIM+1 ellipses, "add a fictional zero row to the end of R",´.
+    % This means that the (fictional) last element of the diagonal is zero 
+    % (i.e. it's index belongs in the list of small elements).
+    % This will trigger the correct solution option later. Don't worry, no panic, read on.
+    if isempty(zeroDiagIndex) &&  numFoci <= DIM
+        zeroDiagIndex = DIM + 1;
+    end
+    
+    % If there are more than one zero on the diagonal, the intersections can not be solved.
+    if length(zeroDiagIndex) > 1
+        return 
+    end
+
+    % If there are no zeros on the diagonal, the problem can be solved
+    if isempty(zeroDiagIndex)
+        xyw = T \ r; % == R \ (Q' * r);
+        intersectionPoints = xyw;
+        return
+    end
+
+    % Now matrix T is made of 3 columns: T = [A B C] and the problem to be solved is
+    % T * [x y w]' = Ax + By + Cw = r  such that  x^2 + y^2 = w^2
+    % Now one of the diagonals of the upper diagonal matrix R is zero,
+    % meaning that the corresponding column of T can be expressed as a linear combination
+    % of the preceeding columns. 
+    % Let's go over all possible positions of the zero diagonal and see what it means.
+    switch zeroDiagIndex
+        case 1
+            % R(1,1) = 0, meaning that column A = 0. 
+            % Solve By + Cw = r and use constraint x^2 = w^2 - y^2
+            yw = T(:,2:3)\r; % Now yw(1) = y, yw(2) = w
+            x = sqrt(yw(2)^2 - yw(1)^2);
+            x = [x  -x];
+            y = [yw(1) yw(1)];
+            w = [yw(2) yw(2)]; % Uncomment if w is needed
+        case 2
+            % R(2,2) = 0, meaning that column B = g * A for some g.
+            % Let's find what g is:
+            g = R(1,2) / R(1,1);
+            % Now we have  Ax + (g * A)y + Cw = r  such that  x^2 + y^2 = w^2
+            % --> A (x + g*y) + Cw = r. 
+            % Denote h=x + g*y, meaning x = h - g*y.
+            % We can now solve h and w without using column B:
+            hw = T(:,[1 3]) \ r; % Now hw(1) = h, hw(2) = w
+            % Now w^2 = x^2 + y^2 = (h - g*y)^2 + y^2
+            % --> (g^2 + 1)*y^2 - 2*g*h*y + (h^2 - w^2) = 0.
+            %  This is a 2nd degree polynomial on y. Let's solve it.
+            y = poly2zeros(g^2 + 1, -2*g*hw(1), hw(1)^2 - hw(2)^2);
+            % Now that we know y, we also know x = -g*y + h
+            x = hw(1) - g*y ;
+            w = [hw(2) hw(2)]; % Uncomment if w is needed
+        case 3
+            % R(3,3) = 0, meaning that column C = g1*A + g2*B for some g1,g2.
+            % Let's find what g is:
+            g = R(1:2,1:2) \ R(1:2,3);
+            % Now A*x + B*y + (g1*A + g2*B)*w = r such that  x^2 + y^2 = w^2
+            % --> A*(x+g1*w) + B*(y+g2*w) = r.
+            % Denote h1 = x+g1*w, h2 = y+g2*w, meaning x = h1 - g1*w, y = h2 - g2*w
+            % Now solve h1 and h2 using columns A and B
+            h = T(:,1:2) \ r;
+            % Now  w^2 = x^2 + y^2 = (h1 - g1*w)^2 + (h2 - g2*w)^2
+            % --> (g1^2 + g2^2 - 1)*w^2 - 2*(g1*h1+g2*h2)*w + (h1^2+h2^2) = 0
+            %  This is a 2nd degree polynomial on w. Let's solve it.
+            w = poly2zeros(g(1)^2 + g(2)^2 - 1, -2*(g(1)*h(1) + g(2)*h(2)), h(1)^2 + h(2)^2);
+            % Now we know w, g and h so we know x and y
+            xy = h - g*w;
+            x = xy(1,:);
+            y = xy(2,:);
+        otherwise
+            error("small diags should be 1 2 or 3 (" + zeroDiagIndex + ")")
+    end
+    intersectionPoints = [x;y;w];
+end
+
+function intersectionPoints = solveSpheroidIntersections3D(fp, diam)
+%   Inputs:
+%       fp - Second focal points as a 3-by-N matrix (the first focal points are always at the origin) 
+%       diam - Diameters of the ellipses along the major semiaxes, N-by-1 vector.
+%   Output:
+%       intersectionPoints - Coordinates of the points where all N ellipses intersect, or an empty vector if they don't intersect.
+    DIM = 3;
+    intersectionPoints = [];
+    numFoci = length(diam); % Number of ellipses (must be DIM or DIM+1)
+
+%  An intersection point (x,y,z) of the spheroids and it's distance w = sqrt(x^2 + y^2 + z^2) from the origin 
+% (i.e. the second focal point) satifies the following conditions:
+%       (x - fp(1,1))^2 + (y - fp(2,1))^2 + (z - fp(3,1))^2 = (diam(1) - w)^2
+%       ...
+%       (x - fp(1,N))^2 + (y - fp(2,N))^2 + (z - fp(3,N))^2 = (diam(N) - w)^2
+%                                           x^2 + y^2 + z^2 = w^2
+% By expanding the squares in N first equations and using the last equation, 
+% we get a linear problem T * [x y z w]' = r  with matrix T and vector r defined as as follows:
+    T = -2 * [fp(1,:)' fp(2,:)' fp(3,:)' diam];
+    r = diam.^2 - fp(1,:)'.^2 - fp(2,:)'.^2 - fp(3,:)'.^2;
+    
+    % QR transform will later help with sorting out corner cases later
+    [~,R] = qr(T);
+    
+    % Are there tiny elements on the diagonal of R? 
+    % If so, it means that some columns of A are linearly dependent.
+    zeroDiagIndex = findSmallIndices(diag(R));
+
+    % If there are fewer than DIM+1 ellipses, "add a fictional zero row to the end of R",´.
+    % This means that the (fictional) last element of the diagonal is zero 
+    % (i.e. it's index belongs in the list of small elements).
+    % This will trigger the correct solution option later. Don't worry, no panic, read on.
+    if isempty(zeroDiagIndex) && numFoci <= DIM
+        zeroDiagIndex = DIM + 1;
+    end
+    
+    % If there are more than one zero on the diagonal, the intersections can not be solved.
+    if length(zeroDiagIndex) > 1
+        return 
+    end
+    
+    % If there are no zeros on the diagonal, the problem can be solved
+    if isempty(zeroDiagIndex)
+        xyzw = T \ r; % == R \ (Q' * r);
+        intersectionPoints = xyzw;
+        return
+    end
+    % Now matrix T is made of 4 columns: T = [A B C D] and the problem to be solved is
+    % T * [x y z w]' = Ax + By + Cz + Dw = r  such that  x^2 + y^2 + z^2 = w^2.
+    % Now one of the diagonals of the upper diagonal matrix R is zero,
+    % meaning that the corresponding column of T can be expressed as a linear combination
+    % of the preceeding columns. 
+    % Let's go over all possible positions of the zero diagonal and see what it means.
+    switch zeroDiagIndex
+        case 1
+            % R(1,1) = 0, meaning that column A = 0. 
+            % Solve By + Cz + Dw = r and use constraint x^2 = w^2 - y^2 - z^2
+            yzw = T(:,2:4)\r; % Now yzw(1) = y, yzw(2) = z, yzw(3) = w
+            x = sqrt(yzw(3)^2 - yzw(1)^2 - yzw(2)^2);
+            x = [x  -x]; % x and -x both are solutions.
+            y = [yzw(1) yzw(1)];
+            z = [yzw(2) yzw(2)];
+            w = [yzw(3) yzw(3)]; % Uncomment if w is needed
+        case 2
+            % R(2,2) = 0, meaning that column B = g * A for some g.
+            % Let's find what g is:
+            g = R(1,2) / R(1,1);
+            % Now we have  Ax + (g * A)y + Cz + Dw = r  such that  x^2 + y^2 + z^2 = w^2
+            % --> A (x + g*y) + Cz + Dw = r. 
+            % Denote h=x + g*y, meaning x = -g*y + h.
+            % We can now solve h and w without using column B:
+            hzw = T(:,[1 3 4]) \ r; % Now hzw(1) = h, hzw(2) = z, hzw(3) = w
+            % Now w^2 = x^2 + y^2 + z^2 = (-g*y + h)^2 + y^2 + z^2
+            % --> (g^2 + 1)*y^2 - 2*g*h*y + (h^2 + z^2 - w^2) = 0.
+            %  This is a 2nd degree polynomial in y. Let's solve it.
+            y = poly2zeros(g^2 + 1, -2*g*hzw(1), hzw(1)^2 + hzw(2)^2 - hzw(3)^2);
+            % Now that we know y, we also know x = -g*y + h
+            x = -g*y + hzw(1);
+            z = [hzw(2) hzw(2)];
+            w = [hzw(3) hzw(3)]; % Uncomment if w is needed
+        case 3
+            % R(3,3) = 0, meaning that column C = g1*A + g2*B for some g1,g2.
+            % Let's find what g is:
+            g = R(1:2,1:2) \ R(1:2,3);
+            % Now A*x + B*y + (g1*A + g2*B)*z + Dw = r such that  x^2 + y^2 + z^2 = w^2
+            % --> A*(x+g1*z) + B*(y+g2*z) + Dw = r.
+            % Denote h1 = x+g1*z, h2 = y+g2*z, meaning x = h1 - g1*z, y = h2 - g2*z
+            % Now solve h1 and h2 using columns A, B and D
+            hw = T(:,[1 2 4]) \ r; % Now h(1:2) are h1,h2 and hw(3) = w
+            % Now w^2 = x^2 + y^2 + z^2 = (h1 - g1*z)^2 + (h2 - g2*z)^2 + z^2
+            % --> (g1^2 + g2^2 + 1)*z^2 - 2*(g1*h1 + g2*h2)*z + (h1^2 + h2^2 - w^2) = 0
+            % This is a 2nd degree polynomial in z. Let's solve it
+            z = poly2zeros(g(1)^2 + g(2)^2 + 1, -2*(g(1)*hw(1) + g(2)*hw(2)), hw(1)^2 + hw(2)^2 - hw(3)^2);
+            xy = hw(1:2) - g*z;
+            x = [xy(1,1) xy(1,2)];
+            y = [xy(2,1) xy(2,2)];
+            w = [hw(3) hw(3)];
+        case 4
+            % R(4,4) = 0, meaning that column D = g1*A + g2*B + g3*C for some g1,g2,g3.
+            % Let's find what g is:
+            g = R(1:3,1:3) \ R(1:3,4);
+            % Now A*x + B*y + C*x + (g1*A + g2*B + g3*C)*w = r such that  x^2 + y^2 +z^2 = w^2
+            % --> A*(x+g1*w) + B*(y+g2*w) + C*(z+g3*w) = r.
+            % Denote h1 = x+g1*w, h2 = y+g2*w, h3 = z+g3*w, 
+            % meaning x = h1 - g1*w, y = h2 - g2*w, z = h3 - g3*w
+            % Now solve h1, h2, h3 using columns A, B and C
+            h = T(:,1:3) \ r;
+
+            % Now w^2 = x^2 + y^2 + z^2 = (h1 - g1*w)^2 + (h2 - g2*w)^2 + (h3 - g3*w)^2
+            % --> (g1^2 + g2^2 + g3^2 - 1)*w^2 - 2*(g1*h1+g2*h2+g3*h3)*w + (h1^2+h2^2+h3^2) = 0
+            %  This is a 2nd degree polynomial on w. Let's solve it.
+            w = poly2zeros(g(1)^2 + g(2)^2 + g(3)^2 - 1, ...
+                           -2*(g(1)*h(1)+g(2)*h(2)+g(3)*h(3)), ...
+                           h(1)^2 + h(2)^2 + h(3)^2);
+            % Now we know w, g and h so we know x, y and z
+            xyz = h - g*w;
+            x = xyz(1,:);
+            y = xyz(2,:);
+            z = xyz(3,:);
+        otherwise
+            error("small diags should be 1 2 or 3 (" + zeroDiagIndex + ")")
+    end
+    intersectionPoints = [x;y;z;w];
+end
+
+function z = poly2zeros(a,b,c)
+% poly2zeros  Finds zeros of polynomial a*x^2 + b*x + c
+    determinant = sqrt(b*b - 4*a*c);
+    z = (-b + [determinant -determinant]) / (2*a);
+end
+
+function smallIndices = findSmallIndices(v)
+% smallIndidices Finds "very small" indices of vector v
+  smallIndices = find(abs(v) <= sqrt(eps(norm(v))));
+end
diff --git a/matlab/threePointDistance.m b/matlab/threePointDistance.m
new file mode 100644
index 0000000..02cd776
--- /dev/null
+++ b/matlab/threePointDistance.m
@@ -0,0 +1,11 @@
+function distance = threePointDistance(point1Array, point2, point3)
+%measuredDistances Calculate distance from each Point1 to Point2 and further on to Point3 
+% Inputs:
+%   point1Array - array of Point1 coordinates, DIM x N array
+%   point2 -  DIM x 1 array
+%   point3 - DIM x 1 array
+% Output:
+%  distances - point1(k) -> point2 -> point3 distance, 1 x N array
+    distance = vecnorm(point1Array - point2, 2, 1);
+    distance = distance + norm(point2 - point3);
+end
\ No newline at end of file
diff --git a/octave/demo2D_2.m b/octave/demo2D_2.m
new file mode 100644
index 0000000..ee9264e
--- /dev/null
+++ b/octave/demo2D_2.m
@@ -0,0 +1,80 @@
+%% This test demonstrates intersections of 2 ellipses
+
+% Coordinate of the shared focal point
+sharedFocalPoint = [-10; -20];
+
+% Coordinates of the other focal points of 2 ellipses.
+numEllipses = 2;
+otherFocalPoints = zeros(2, numEllipses);
+
+%% (Case 1: The major axes are collinear and upright)
+otherFocalPoints(:,1) = [0; 20] + sharedFocalPoint;
+otherFocalPoints(:,2) = [0; 10] + sharedFocalPoint;
+
+% Select point p at which the ellipses intersect and determine the diameters accordingly.
+% There will be 2 points of intersection, one of which will be p.
+p = [4; 3] + sharedFocalPoint;
+diameters = threePointDistance(otherFocalPoints, p, sharedFocalPoint)';
+
+% Solve intersections and an estimate on the model error.
+% The bigger the error, the further away the estimated intersections are from the real ones. 
+% The model error is not exacly the actual distance between 
+% the estimate and the correct intersection so it should be taken with a pinch of salt.
+[intersections, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters);
+
+assert(min(vecnorm(intersections - p)) < 1e-6, "Intersection point is not correct")
+assert(max(relativeModelError) < 1e-6, "Model error in the solver is too large")
+
+% Calculate ellipse coordinates and plot the coordinates and the intersection points.
+plotEllipses(sharedFocalPoint, otherFocalPoints, diameters, intersections, "2 ellipses, Case 1");
+
+disp("Case 1: intersection coordinates as columns: ");
+disp(intersections);
+
+%% (Case 2: The major axes are collinear and tilted)
+otherFocalPoints(:,1) = [-10; 20] + sharedFocalPoint;
+otherFocalPoints(:,2) = [-5; 10] + sharedFocalPoint;
+
+% Select point p at which the ellipses intersect and determine the diameters accordingly.
+% There will be 2 points of intersection, one of which will be p.
+p = [4; 3] + sharedFocalPoint;
+diameters = threePointDistance(otherFocalPoints, p, sharedFocalPoint)';
+
+% Solve intersections and an estimate on the model error.
+% The bigger the error, the further away the estimated intersections are from the real ones. 
+% The model error is not exacly the actual distance between 
+% the estimate and the correct intersection so it should be taken with a pinch of salt.
+[intersections, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters);
+
+assert(min(vecnorm(intersections - p)) < 1e-6, "Intersection point is not correct")
+assert(max(relativeModelError) < 1e-6, "Model error in the solver is too large")
+
+% Calculate ellipse coordinates and plot the coordinates and the intersection points.
+plotEllipses(sharedFocalPoint, otherFocalPoints, diameters, intersections, "2 ellipses, Case 2");
+
+disp("Case 2: intersection coordinates as columns:");
+disp(intersections);
+
+%% (Case 3: General case with 2 ellipses)
+otherFocalPoints(:,1) = [10; 4] + sharedFocalPoint;
+otherFocalPoints(:,2) = [0; 10] + sharedFocalPoint;
+
+% Select point p at which the ellipses intersect and determine the diameters accordingly.
+% There will be 2 points of intersection, one of which will be p.
+p = [4; 3] + sharedFocalPoint;
+diameters = threePointDistance(otherFocalPoints, p, sharedFocalPoint)';
+
+% Solve intersections and an estimate on the model error.
+% The bigger the error, the further away the estimated intersections are from the real ones. 
+% The model error is not exacly the actual distance between 
+% the estimate and the correct intersection so it should be taken with a pinch of salt.
+[intersections, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters);
+
+assert(min(vecnorm(intersections - p)) < 1e-6, "Intersection point is not correct")
+assert(max(relativeModelError) < 1e-6, "Model error in the solver is too large")
+
+% Calculate ellipse coordinates and plot the coordinates and the intersection points.
+plotEllipses(sharedFocalPoint, otherFocalPoints, diameters, intersections, "2 ellipses, Case 3");
+
+disp("Case 3: intersection coordinates as columns:");
+disp(intersections);
\ No newline at end of file
diff --git a/octave/demo2D_3.m b/octave/demo2D_3.m
new file mode 100644
index 0000000..9857215
--- /dev/null
+++ b/octave/demo2D_3.m
@@ -0,0 +1,112 @@
+%% This test demonstrates intersections of 3 ellipses
+
+% Coordinate of the shared focal point
+sharedFocalPoint = [-10; -20];
+
+% Coordinates of the other focal points of 2 ellipses.
+numEllipses = 3;
+otherFocalPoints = zeros(2, numEllipses);
+
+%% (Case 1: The major axes are collinear and upright)
+otherFocalPoints(:,1) = [0; 20] + sharedFocalPoint;
+otherFocalPoints(:,2) = [0; 10] + sharedFocalPoint;
+otherFocalPoints(:,3) = [0; -16] + sharedFocalPoint;
+
+% Select point p at which the ellipses intersect and determine the diameters accordingly.
+% There will be 2 points of intersection, one of which will be p.
+p = [4; 3] + sharedFocalPoint;
+diameters = threePointDistance(otherFocalPoints, p, sharedFocalPoint)';
+
+% Solve intersections and an estimate on the model error.
+% The bigger the error, the further away the estimated intersections are from the real ones.
+% The model error is not exacly the actual distance between
+% the estimate and the correct intersection so it should be taken with a pinch of salt.
+[intersections, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters);
+
+assert(min(vecnorm(intersections - p)) < 1e-6, "Intersection point is not correct")
+assert(max(relativeModelError) < 1e-6, "Model error in the solver is too large")
+
+% Calculate ellipse coordinates and plot the coordinates and the intersection points.
+plotEllipses(sharedFocalPoint, otherFocalPoints, diameters, intersections, "3 ellipses, Case 1");
+
+disp("Case 1: intersection coordinates as columns: ");
+disp(intersections);
+
+%% (Case 2: The major axes are collinear and tilted)
+otherFocalPoints(:,1) = [-10; 20] + sharedFocalPoint;
+otherFocalPoints(:,2) = [-5; 10] + sharedFocalPoint;
+otherFocalPoints(:,3) = [-8; 16] + sharedFocalPoint;
+
+% Select point p at which the ellipses intersect and determine the diameters accordingly.
+% There will be 2 points of intersection, one of which will be p.
+p = [4; 3] + sharedFocalPoint;
+diameters = threePointDistance(otherFocalPoints, p, sharedFocalPoint)';
+
+% Solve intersections and an estimate on the model error.
+% The bigger the error, the further away the estimated intersections are from the real ones.
+% The model error is not exacly the actual distance between
+% the estimate and the correct intersection so it should be taken with a pinch of salt.
+[intersections, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters);
+
+assert(min(vecnorm(intersections - p)) < 1e-6, "Intersection point is not correct")
+assert(max(relativeModelError) < 1e-6, "Model error in the solver is too large")
+
+% Calculate ellipse coordinates and plot the coordinates and the intersection points.
+plotEllipses(sharedFocalPoint, otherFocalPoints, diameters, intersections, "3 ellipses, Case 2");
+
+disp("Case 2: intersection coordinates as columns:");
+disp(intersections);
+
+%% (Case 3: General case with 3 ellipses)
+otherFocalPoints(:,1) = [10; 4] + sharedFocalPoint;
+otherFocalPoints(:,2) = [0; 10] + sharedFocalPoint;
+otherFocalPoints(:,3) = [10; 0] + sharedFocalPoint;
+
+% Select point p at which the ellipses intersect and determine the diameters accordingly.
+% There will be 2 points of intersection, one of which will be p.
+p = [4; 3] + sharedFocalPoint;
+diameters = threePointDistance(otherFocalPoints, p, sharedFocalPoint)';
+
+% Solve intersections and an estimate on the model error.
+% The bigger the error, the further away the estimated intersections are from the real ones.
+% The model error is not exacly the actual distance between
+% the estimate and the correct intersection so it should be taken with a pinch of salt.
+[intersections, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters);
+
+assert(min(vecnorm(intersections - p)) < 1e-6, "Intersection point is not correct")
+assert(max(relativeModelError) < 1e-6, "Model error in the solver is too large")
+
+% Calculate ellipse coordinates and plot the coordinates and the intersection points.
+plotEllipses(sharedFocalPoint, otherFocalPoints, diameters, intersections, "3 ellipses, Case 3");
+
+disp("Case 3: (One intersection only) Intersection coordinates as columns:");
+disp(intersections);
+
+%% (Case 4: Non-solvable case with 3 ellipses)
+% This is an example of the solver giving a wrong answer and returning a high
+% model error when the requirement that the ellipses intersect at one point does not hold.
+otherFocalPoints(:,1) = [10; 4] + sharedFocalPoint;
+otherFocalPoints(:,2) = [0; 10] + sharedFocalPoint;
+otherFocalPoints(:,3) = [10; 0] + sharedFocalPoint;
+
+% Select point p at which the ellipses intersect and determine the diameters accordingly.
+% There will be 2 points of intersection, one of which will be p.
+p = [4; 3] + sharedFocalPoint;
+diameters = threePointDistance(otherFocalPoints, p, sharedFocalPoint)';
+
+% The 3rd ellipse does not intersect at the same point where the 1st and 2nd ellipses intersect
+pp = p + [1; 1];
+diameters(1) = threePointDistance(otherFocalPoints(:,1), pp, sharedFocalPoint)';
+
+% Solve intersections and an estimate on the model error.
+% The bigger the error, the further away the estimated intersections are from the real ones.
+% The model error is not exacly the actual distance between
+% the estimate and the correct intersection so it should be taken with a pinch of salt.
+[intersections, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters);
+
+% Calculate ellipse coordinates and plot the coordinates and the intersection points.
+plotEllipses(sharedFocalPoint, otherFocalPoints, diameters, intersections, ...
+             ["3 ellipses intersecting at different points, rel.model error=" num2str(relativeModelError)]);
+
+disp(["Case 4: intersection with model error " num2str(relativeModelError)]);
+disp(intersections);
diff --git a/octave/demo3D_3.m b/octave/demo3D_3.m
new file mode 100644
index 0000000..d2525f1
--- /dev/null
+++ b/octave/demo3D_3.m
@@ -0,0 +1,107 @@
+%% This test demonstrates intersections of 3 spheroids
+
+% Coordinate of the shared focal point
+sharedFocalPoint = [-10; -20; 30];
+
+% Coordinates of the other focal points of 2 spheroids.
+numSpheroids = 3;
+otherFocalPoints = zeros(3, numSpheroids);
+
+%% (Case 1: The major axes are on the same yz-plane where x is constant)
+otherFocalPoints(:,1) = [0; 20; 1] + sharedFocalPoint;
+otherFocalPoints(:,2) = [0; 10; -2] + sharedFocalPoint;
+otherFocalPoints(:,3) = [0; -16; 3] + sharedFocalPoint;
+
+% Select point p at which the spheroids intersect and determine the diameters accordingly.
+% There will be 2 points of intersection, one of which will be p.
+p = [4; 3; 11] + sharedFocalPoint;
+diameters = threePointDistance(otherFocalPoints, p, sharedFocalPoint)';
+
+% Solve intersections and an estimate on the model error.
+% The bigger the error, the further away the estimated intersections are from the real ones. 
+% The model error is not exacly the actual distance between 
+% the estimate and the correct intersection so it should be taken with a pinch of salt.
+[intersections, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters);
+
+assert(min(vecnorm(intersections - p)) < 1e-6, "Intersection point is not correct")
+assert(max(relativeModelError) < 1e-6, "Model error in the solver is too large")
+
+% Calculate spheroid coordinates and plot the coordinates and the intersection points.
+plotSpheroids(sharedFocalPoint, otherFocalPoints, diameters, intersections, "3 spheroids, Case 1");
+
+disp("Case 1: intersection coordinates as columns: ");
+disp(intersections);
+
+%% (Case 2: The major axes are on the same tilted plane (x and y coordinates are linearly dependent))
+otherFocalPoints(:,1) = [10; 20; 1] + sharedFocalPoint;
+otherFocalPoints(:,2) = [5; 10; -2] + sharedFocalPoint;
+otherFocalPoints(:,3) = [-8; -16; 3] + sharedFocalPoint;
+
+% Select point p at which the spheroids intersect and determine the diameters accordingly.
+% There will be 2 points of intersection, one of which will be p.
+p = [4; 3; 11] + sharedFocalPoint;
+diameters = threePointDistance(otherFocalPoints, p, sharedFocalPoint)';
+
+% Solve intersections and an estimate on the model error.
+% The bigger the error, the further away the estimated intersections are from the real ones. 
+% The model error is not exacly the actual distance between 
+% the estimate and the correct intersection so it should be taken with a pinch of salt.
+[intersections, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters);
+
+assert(min(vecnorm(intersections - p)) < 1e-6, "Intersection point is not correct")
+assert(max(relativeModelError) < 1e-6, "Model error in the solver is too large")
+
+% Calculate spheroid coordinates and plot the coordinates and the intersection points.
+plotSpheroids(sharedFocalPoint, otherFocalPoints, diameters, intersections, "3 spheroids, Case 2");
+
+disp("Case 2: intersection coordinates as columns:");
+disp(intersections);
+
+%% (Case 3: z coordinates of the major axes are linearly dependent on x- and y-coordinates.)
+otherFocalPoints(:,1) =[10; 1; 10+1] + sharedFocalPoint;
+otherFocalPoints(:,2) = [5; -2; 5-2] + sharedFocalPoint;
+otherFocalPoints(:,3) = [-8; 3; -8+3] + sharedFocalPoint;
+
+% Select point p at which the spheroids intersect and determine the diameters accordingly.
+% There will be 2 points of intersection, one of which will be p.
+p = [4; 3; 11] + sharedFocalPoint;
+diameters = threePointDistance(otherFocalPoints, p, sharedFocalPoint)';
+
+% Solve intersections and an estimate on the model error.
+% The bigger the error, the further away the estimated intersections are from the real ones. 
+% The model error is not exacly the actual distance between 
+% the estimate and the correct intersection so it should be taken with a pinch of salt.
+[intersections, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters);
+
+assert(min(vecnorm(intersections - p)) < 1e-6, "Intersection point is not correct")
+assert(max(relativeModelError) < 1e-6, "Model error in the solver is too large")
+
+% Calculate spheroid coordinates and plot the coordinates and the intersection points.
+plotSpheroids(sharedFocalPoint, otherFocalPoints, diameters, intersections, "3 spheroids, Case 3");
+
+disp("Case 3: intersection coordinates as columns:");
+disp(intersections);
+%% (Case 4: General case with 3 spheroids)
+otherFocalPoints(:,1) = [-2; -5; 1] + sharedFocalPoint;
+otherFocalPoints(:,2) = [0; 10; 2] + sharedFocalPoint;
+otherFocalPoints(:,3) = [10; 0; 3] + sharedFocalPoint;
+
+% Select point p at which the spheroids intersect and determine the diameters accordingly.
+% There will be 2 points of intersection, one of which will be p.
+p = [4; 3; 11] + sharedFocalPoint;
+diameters = threePointDistance(otherFocalPoints, p, sharedFocalPoint)';
+
+% Solve intersections and an estimate on the model error.
+% The bigger the error, the further away the estimated intersections are from the real ones. 
+% The model error is not exacly the actual distance between 
+% the estimate and the correct intersection so it should be taken with a pinch of salt.
+[intersections, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters);
+
+assert(min(vecnorm(intersections - p)) < 1e-6, "Intersection point is not correct")
+assert(max(relativeModelError) < 1e-6, "Model error in the solver is too large")
+
+% Calculate spheroid coordinates and plot the coordinates and the intersection points.
+plotSpheroids(sharedFocalPoint, otherFocalPoints, diameters, intersections, "3 spheroids, Case 4");
+
+disp("Case 4: intersection coordinates as columns:");
+disp(intersections);
\ No newline at end of file
diff --git a/octave/demo3D_4.m b/octave/demo3D_4.m
new file mode 100644
index 0000000..e58f727
--- /dev/null
+++ b/octave/demo3D_4.m
@@ -0,0 +1,111 @@
+%% This test demonstrates intersections of 4 spheroids
+
+% Coordinate of the shared focal point
+sharedFocalPoint = [-10; -20; 30];
+
+% Coordinates of the other focal points of 2 spheroids.
+numSpheroids = 4;
+otherFocalPoints = zeros(3, numSpheroids);
+
+%% (Case 1: The major axes are on the same yz-plane where x is constant)
+otherFocalPoints(:,1) = [0; 20; 1] + sharedFocalPoint;
+otherFocalPoints(:,2) = [0; 10; -2] + sharedFocalPoint;
+otherFocalPoints(:,3) = [0; -16; 3] + sharedFocalPoint;
+otherFocalPoints(:,4) = [0; -14; -4] + sharedFocalPoint;
+
+% Select point p at which the spheroids intersect and determine the diameters accordingly.
+% There will be 2 points of intersection, one of which will be p.
+p = [4; 3; 11] + sharedFocalPoint;
+diameters = threePointDistance(otherFocalPoints, p, sharedFocalPoint)';
+
+% Solve intersections and an estimate on the model error.
+% The bigger the error, the further away the estimated intersections are from the real ones. 
+% The model error is not exacly the actual distance between 
+% the estimate and the correct intersection so it should be taken with a pinch of salt.
+[intersections, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters);
+
+assert(min(vecnorm(intersections - p)) < 1e-6, "Intersection point is not correct")
+assert(max(relativeModelError) < 1e-6, "Model error in the solver is too large")
+
+% Calculate spheroid coordinates and plot the coordinates and the intersection points.
+plotSpheroids(sharedFocalPoint, otherFocalPoints, diameters, intersections, "4 spheroids, Case 1");
+
+disp("Case 1: intersection coordinates as columns: ");
+disp(intersections);
+
+%% (Case 2: The major axes are on the same tilted plane (x and y coordinates are linearly dependent))
+otherFocalPoints(:,1) = [10; 20; 1] + sharedFocalPoint;
+otherFocalPoints(:,2) = [5; 10; -2] + sharedFocalPoint;
+otherFocalPoints(:,3) = [-8; -16; 3] + sharedFocalPoint;
+otherFocalPoints(:,4) = [-7; -14; -4] + sharedFocalPoint;
+
+% Select point p at which the spheroids intersect and determine the diameters accordingly.
+% There will be 2 points of intersection, one of which will be p.
+p = [4; 3; 11] + sharedFocalPoint;
+diameters = threePointDistance(otherFocalPoints, p, sharedFocalPoint)';
+
+% Solve intersections and an estimate on the model error.
+% The bigger the error, the further away the estimated intersections are from the real ones. 
+% The model error is not exacly the actual distance between 
+% the estimate and the correct intersection so it should be taken with a pinch of salt.
+[intersections, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters);
+
+assert(min(vecnorm(intersections - p)) < 1e-6, "Intersection point is not correct")
+assert(max(relativeModelError) < 1e-6, "Model error in the solver is too large")
+
+% Calculate spheroid coordinates and plot the coordinates and the intersection points.
+plotSpheroids(sharedFocalPoint, otherFocalPoints, diameters, intersections, "4 spheroids, Case 2");
+
+disp("Case 2: intersection coordinates as columns:");
+disp(intersections);
+
+%% (Case 3: z coordinates of the major axes are linearly dependent on x- and y-coordinates.)
+otherFocalPoints(:,1) =[10; 1; 10+1] + sharedFocalPoint;
+otherFocalPoints(:,2) = [5; -2; 5-2] + sharedFocalPoint;
+otherFocalPoints(:,3) = [-8; 3; -8+3] + sharedFocalPoint;
+otherFocalPoints(:,4) = [-7; -4; -7-4] + sharedFocalPoint;
+
+% Select point p at which the spheroids intersect and determine the diameters accordingly.
+% There will be 2 points of intersection, one of which will be p.
+p = [4; 3; 11] + sharedFocalPoint;
+diameters = threePointDistance(otherFocalPoints, p, sharedFocalPoint)';
+
+% Solve intersections and an estimate on the model error.
+% The bigger the error, the further away the estimated intersections are from the real ones. 
+% The model error is not exacly the actual distance between 
+% the estimate and the correct intersection so it should be taken with a pinch of salt.
+[intersections, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters);
+
+assert(min(vecnorm(intersections - p)) < 1e-6, "Intersection point is not correct")
+assert(max(relativeModelError) < 1e-6, "Model error in the solver is too large")
+
+% Calculate spheroid coordinates and plot the coordinates and the intersection points.
+plotSpheroids(sharedFocalPoint, otherFocalPoints, diameters, intersections, "4 spheroids, Case 3");
+
+disp("Case 3: intersection coordinates as columns:");
+disp(intersections);
+%% (Case 4: General case with 3 spheroids)
+otherFocalPoints(:,1) = [-2; -5; 1] + sharedFocalPoint;
+otherFocalPoints(:,2) = [0; 10; 2] + sharedFocalPoint;
+otherFocalPoints(:,3) = [10; 0; 3] + sharedFocalPoint;
+otherFocalPoints(:,4) = [-8; -3; 4] + sharedFocalPoint;
+
+% Select point p at which the spheroids intersect and determine the diameters accordingly.
+% There will be 2 points of intersection, one of which will be p.
+p = [4; 3; 11] + sharedFocalPoint;
+diameters = threePointDistance(otherFocalPoints, p, sharedFocalPoint)';
+
+% Solve intersections and an estimate on the model error.
+% The bigger the error, the further away the estimated intersections are from the real ones. 
+% The model error is not exacly the actual distance between 
+% the estimate and the correct intersection so it should be taken with a pinch of salt.
+[intersections, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters);
+
+assert(min(vecnorm(intersections - p)) < 1e-6, "Intersection point is not correct")
+assert(max(relativeModelError) < 1e-6, "Model error in the solver is too large")
+
+% Calculate spheroid coordinates and plot the coordinates and the intersection points.
+plotSpheroids(sharedFocalPoint, otherFocalPoints, diameters, intersections, "4 spheroids, Case 4");
+
+disp("Case 4: intersection coordinates as columns:");
+disp(intersections);
\ No newline at end of file
diff --git a/octave/plotEllipses.m b/octave/plotEllipses.m
new file mode 100644
index 0000000..1812cff
--- /dev/null
+++ b/octave/plotEllipses.m
@@ -0,0 +1,82 @@
+function output = plotEllipses(sharedFocalPoint, otherFocalPoints, diameters, intersections, titleText)
+%plotEllipses Plot ellipses with the given focal points and diameters. Also show the points of intersections.
+% Inputs:
+%   sharedFocalPoint - Coordinate of shared focal points of every ellipse (2x1)
+%   otherFocalPoints - Coordinates of N other focal points (2xN)
+%   diameters - Diameters of the ellipses along the major axis (Nx1)
+%   intersections - Coordinates of M intersections (2xM)
+
+    if nargin < 5
+        titleText = "";
+    end
+
+    numEllipses = size(otherFocalPoints, 2);
+    numPoints = 100;
+    ellipseCoord = zeros(2,numPoints, numEllipses);
+    theta = linspace(0, 2*pi, numPoints);
+    diffCoord = otherFocalPoints - sharedFocalPoint;
+    for ii = 1:numEllipses
+        a = diameters(ii)/2;  % Major axis
+        b = sqrt(a*a - norm(diffCoord(:,ii)/2)^2); % Minor axis
+        phase = diffCoord(1,ii) + 1i*diffCoord(2,ii);
+        phase = phase / norm(phase); % Rotation of the ellipse
+        center = (otherFocalPoints(:,ii) + sharedFocalPoint) / 2; % Center of the ellipse
+        for jj = 1:numPoints
+            coord = (a * cos(theta(jj))) + 1i*(b * sin(theta(jj))); % axis aligned
+            coord = coord * phase;  % rotate
+            coord = coord + (center(1) + 1i*center(2));
+            ellipseCoord(:,jj,ii) = [real(coord); imag(coord)];
+        end
+    end
+    rgb=['g','r','b'];
+    limX = [floor(min(ellipseCoord(1,:))) ceil(max(ellipseCoord(1,:)))];
+    limY = [floor(min(ellipseCoord(2,:))) ceil(max(ellipseCoord(2,:)))];
+    legendArray = {};
+
+    figure;
+    hold on
+    % Plot the shared focal point
+    plot(sharedFocalPoint(1),sharedFocalPoint(2),'Marker','o','MarkerSize',10,'Color','k');
+    legendArray{end+1} = "Shared focal point";
+
+    % Plot non-shared focal points
+    for ii=1:numEllipses
+        plot(otherFocalPoints(1,ii),otherFocalPoints(2,ii),'Marker','*','MarkerSize',10,'Color',rgb(mod(ii,3)+1));
+        legendArray{end+1} = ["focal point #"  num2str(ii)];
+    end
+
+    % Plot ellipses
+    for ii=1:numEllipses
+        plot(ellipseCoord(1,:,ii), ellipseCoord(2,:,ii),'-','Color',rgb(mod(ii,3)+1));
+        legendArray{end+1} = ["Ellipse #"  num2str(ii)];
+    end
+
+
+    % Plot the intersections
+    for jj=1:size(intersections,2)
+        if mod(jj,2) == 0
+            mkr = 'pentagram';
+        else
+            mkr = 'hexagram';
+        end
+        plot(intersections(1,jj), intersections(2,jj),'Marker',mkr,'MarkerSize',12,'Color','k','LineStyle',':');
+        legendArray{end+1} = ["Intersection #" num2str(jj)];
+    end
+
+    for jj=1:size(intersections,2)
+        for ii=1:numEllipses
+            plot([intersections(1,jj) otherFocalPoints(1,ii)], [intersections(2,jj) otherFocalPoints(2,ii)],':','Color',rgb(mod(ii,3)+1))
+            plot([intersections(1,jj) sharedFocalPoint(1)], [intersections(2,jj) sharedFocalPoint(2)],':','Color',rgb(mod(ii,3)+1))
+        end
+    end
+
+    grid on
+    axis equal
+    ax = gca;
+    set (ax, "XTick", limX(1):limX(2));
+    set (ax, "YTick", limY(1):limY(2));
+    title(titleText)
+    legend(legendArray);
+
+    output = gca;
+end
diff --git a/octave/plotSpheroids.m b/octave/plotSpheroids.m
new file mode 100644
index 0000000..7bc7ed1
--- /dev/null
+++ b/octave/plotSpheroids.m
@@ -0,0 +1,96 @@
+function output = plotSpheroids(sharedFocalPoint, otherFocalPoints, diameters, intersections, titleText)
+%plotEllipses Plot spheroids with the given focal points and diameters. Also show the points of intersections.
+% Inputs:
+%   sharedFocalPoint - Coordinate of shared focal points of every ellipse (2x1)
+%   otherFocalPoints - Coordinates of N other focal points (2xN)
+%   diameters - Diameters of the ellipses along the major axis (Nx1)
+%   intersections - Coordinates of M intersections (2xM)
+
+    if nargin < 5
+        titleText = "";
+    end
+
+    numPoints = 50;
+    numSpheroids = size(otherFocalPoints, 2);
+    spheroidCoord = zeros(3,numPoints, numSpheroids);
+    theta = linspace(0, 2*pi, numPoints);
+    diffCoord = otherFocalPoints - sharedFocalPoint;
+
+    figure;
+    hold on
+    grid on
+    axis equal
+    legendArray = {};
+
+    for ii = 1:numSpheroids
+        a = diameters(ii)/2;  % Major axis
+        b = sqrt(a*a - norm(diffCoord(:,ii)/2)^2); % Minor axis
+        helperEllipse = zeros(3, numPoints*numPoints);
+        helperIndex = 0;
+        for jj = 1:numPoints
+            coord = [a * cos(theta(jj)) , b * sin(theta(jj))];
+            for kk = 1:numPoints
+                helperIndex = helperIndex + 1;
+                helperEllipse(:,helperIndex) = [coord(1); cos(theta(kk))*coord(2); sin(theta(kk))*coord(2)];
+            end
+        end
+
+        center = (otherFocalPoints(:,ii) + sharedFocalPoint) / 2; % Center of the spheroid
+
+        % Cartesian x axis will be mapped into this vector,
+        % which is the direction of the major axis.
+        Vaxis = diffCoord(:,ii);
+        Vaxis = Vaxis / norm(Vaxis);
+        % Decide another axis orthogonal the Xaxis
+        [~,minIndex] = min(abs(axis));
+        switch minIndex
+            case 1
+                Uaxis = [1;0;0];
+            case 2
+                Uaxis = [0;1;0];
+            case 3
+                Uaxis = [0;0;1];
+        end
+        Waxis = Uaxis - (Uaxis' * Vaxis) * Vaxis;
+        Waxis = Waxis / norm(Waxis);
+        Uaxis=cross(Waxis, Vaxis);
+        % This matrix rotates the helper ellipse into place
+        Rot = [Vaxis Waxis Uaxis];
+
+        spheroid = Rot * helperEllipse + center;
+
+        plot3(spheroid(1,:), spheroid(2,:), spheroid(3,:),'.','MarkerSize',10);
+        legendArray{end+1} = ["Spheroid " num2str(ii)];
+        % Plot Rx->Tx pairs
+        plot3([sharedFocalPoint(1), otherFocalPoints(1,ii)], ...
+            [sharedFocalPoint(2), otherFocalPoints(2,ii)], ...
+            [sharedFocalPoint(3), otherFocalPoints(3,ii)],'x-','LineWidth',3);
+        legendArray{end+1} = ["Foci axis " num2str(ii)];
+    end
+
+    % Plot the intersections
+    numIntersections = size(intersections,2);
+    if numIntersections == 1
+            plot3(intersections(1), intersections(2), intersections(3), ...
+                'ro', 'MarkerSize',16,'MarkerFaceColor','auto','LineWidth',2);
+            legendArray{end+1} = "Intersection";
+    elseif numIntersections > 1
+        for jj=1:size(intersections,2)
+            if mod(jj,2) == 0
+                mkr = 'pentagram';
+            else
+                mkr = 'hexagram';
+            end
+            plot3(intersections(1,jj), intersections(2,jj), intersections(3,jj), ...
+                'Color','r','Marker',mkr, 'MarkerSize',16,'MarkerFaceColor','auto','LineWidth',2);
+            legendArray{end+1} = ["Intersection #"  num2str(jj)];
+        end
+    end
+
+    xlabel("XX")
+    ylabel("YY")
+    zlabel("ZZ")
+    legend(legendArray)
+    title(titleText)
+    output = gca;
+end
diff --git a/octave/solveEllipseIntersections.m b/octave/solveEllipseIntersections.m
new file mode 100644
index 0000000..4c955f5
--- /dev/null
+++ b/octave/solveEllipseIntersections.m
@@ -0,0 +1,335 @@
+function [intersectionPoints, relativeModelError] = solveEllipseIntersections(sharedFocalPoint, otherFocalPoints, diameters)
+%solveIntersections Calculates intersections of ellipses and spheroids with a shared focal point.
+%   Calculates intersections of ellipses (2D) or spheroids (3D). All ellipses are assumed to have the first focal point  
+%   at the same position and the second focal point at the given coordinates.
+%   The algorithm for 3D spheroids is described below. The 2D case is similar but simpler.
+%
+%   The diameters along the major semiaxes determine the sum of distances from the first focal point to 
+%   any point on the locus of the ellipse or spheroid and further on to the second focal point.
+%   Without loss of generality we assume that the first focal point of each spheroid is at the origin 
+%   and fp = otherFocalPoints - sharedFocalPoint is the set of the second focal points.
+%   Now assume that [x y z] is an intersection of the spheroids. Its distance from the origin is denoted by w.
+%   Now the intersection point must satisfy the following N+1 conditions where N is the number of intersecting spheroids:
+%
+%                                           x^2 + y^2 + z^2 = w^2
+%       (x - fp(1,1))^2 + (y - fp(2,1))^2 + (z - fp(3,1))^2 = (diam(1) - w)^2
+%       ...
+%       (x - fp(1,N))^2 + (y - fp(2,N))^2 + (z - fp(3,N))^2 = (diam(N) - w)^2
+%
+% By expanding the squares in N last equations and subtracting the first equation
+% from the expanded squares, we get an N-by-4 linear system  T*[x y z w]' = r as follows:
+%
+%   [-2 fp(1,1) -2 fp(2,1) -2 fp(3,1) diam(1)] [x] = [diam(1)^2 - fp(1,1)^2 - fp(2,1)^2 - fp(3,1)^2]
+%   [                   ...                  ] [y] = [      ....                       ]
+%   [                   ...                  ] [z] = [      ....                       ]
+%   [-2 fp(1,N) -2 fp(2,N) -2 fp(3,N) diam(N)] [w] = [diam(N)^2 - fp(1,N)^2 - fp(2,N)^2 - fp(3,N)^2]
+%
+%   So T = -2 * [fp(1,:)' fp(2,:)' fp(3,:)' diam] and r =  diam.^2 - fp(1,:)'.^2 - fp(2,:)'.^2 - fp(3,:)'.^2;
+%
+%   Notice that a row of T may be linearly dependent on the other rows. Those cases can be dealt with
+%   by writing the dependent row as a linear combination of the other rows and using the fact that
+%   x^2 + y^2 + z^2 = w^2 to solve the coefficients of the linear combination. This will be explained in
+%   comments in the code below.
+%
+%   Inputs:
+%       sharedFocalPoint - The focal point which is shared by all ellipses / spheroids, DIM-by-1 vector
+%       otherFocalPoints - The other (non-shared) focal points as DIM-by-N matrix.
+%          In 2D case, N should be 2 or 3. In 3D case, N should be 3 or 4.
+%       diameters - Diameters of the ellipses along the major semiaxes as an N-by-1 vector.
+%   Output:
+%       intersectionPoints - Coordinates of the points where all N ellipses intersect. 
+%                            A 2-by-M or 3-by-M matrix where M is 1 of the intersection point is unique and 2 if there are two candidates.
+%                            If the matrix is empty, there is no common intersection point.
+%                            Note that if there is noise or other inaccuracy in the set of diameters, the all ellipses/spheroids
+%                            may not actually intersect at the given points. In this case distanceError will be large
+%                            so the inaccuracy can be detected.
+%       relativeModelError - Relative difference between the actual squared distance of the calculated intersection point(s)
+%                            from the shared focal point (i.e. x^2+y^2+z^2) and the squared distance predicted by
+%                            the model (i.e. w^2). So relativeModelError = sqrt(abs(x^2+y^2+z^2 - w^2) / x^2+y^2+z^2).
+%                            The larger the value, the less reliable the result is. 
+%                            The value is very near zero if all ellipses/spheroids intersect at the same points.
+
+assert(nargin == 3)
+% 2D or 3D mode?
+
+dimension = size(otherFocalPoints, 1);
+numFoci = size(otherFocalPoints, 2);
+
+assert(numFoci == length(diameters), "The number of focal points and the number of diameters must be the same.");
+assert(size(diameters, 2) == 1, "Diameters must be an N-by-1 vector")
+assert(size(sharedFocalPoint,1) == dimension)
+assert(size(sharedFocalPoint,2) == 1)
+
+% Normalize the coordinates so that the shared focal point is at the origin.
+focalPoints = otherFocalPoints - sharedFocalPoint;
+switch dimension
+    case 2
+        assert(numFoci == 2 || numFoci == 3, "In 2D case the number of points should be 2 or 3")
+         intersectionPoints = solveEllipseIntersections2D(focalPoints, diameters);
+    case 3
+        assert(numFoci == 3 || numFoci == 4, "In 3D case the number of points should be 3 or 4")
+        intersectionPoints = solveSpheroidIntersections3D(focalPoints, diameters);
+    otherwise
+        error("The dimension of focal points must be 2 or 3")
+end
+
+
+if isempty(intersectionPoints)
+    relativeModelError = Inf;
+    return
+end
+
+% Sometimes the imaginary components can be non-zero due to numerical errors in the solver.
+if ~isreal(intersectionPoints)
+    % The points are conjugates, so only one real part is enough.
+    intersectionPoints = real(intersectionPoints(:,1));
+end
+if nargout > 1
+    w = intersectionPoints(dimension+1,:);
+    xyz = intersectionPoints(1:dimension,:);
+    relativeModelError =  sqrt(abs(vecnorm(xyz).^2 - w.^2)) ./ (vecnorm(xyz) + realmin);
+end
+
+% Undo the normalization
+intersectionPoints = intersectionPoints(1:dimension,:) + sharedFocalPoint;
+end
+
+
+function intersectionPoints = solveEllipseIntersections2D(fp, diam)
+%   Inputs:
+%       fp - Second focal points as a 2-by-N matrix (the first focal points are always at the origin) 
+%       diam - Diameters of the ellipses along the major semiaxes, N-by-1 vector.
+%   Output:
+%       intersectionPoints - Coordinates of the points where all N ellipses intersect, or an empty vector if they don't intersect.
+%
+    DIM = 2;
+    intersectionPoints = [];
+    numFoci = length(diam); % Number of ellipses (must be DIM or DIM+1)
+
+%  An intersection point (x,y) of the spheroids and it's distance w = sqrt(x^2 + y^2) from the origin 
+% (i.e. the second focal point) satifies the following conditions:
+%       (x - fp(1,1))^2 + (y - fp(2,1))^2 = (diam(1) - w)^2
+%       ...
+%       (x - fp(1,N))^2 + (y - fp(2,N))^2 = (diam(N) - w)^2
+%                               x^2 + y^2 = w^2
+% By expanding the squares in N first equations and using the last equation, 
+% we get a linear problem T * [x y w]' = r  with matrix T and vector r defined as as follows:
+    T = -2 * [fp(1,:)' fp(2,:)' diam];
+    r = diam.^2 - fp(1,:)'.^2 - fp(2,:)'.^2;
+    
+    % QR transform will later help with sorting out corner cases later
+    [~,R] = qr(T);
+    
+    % Are there tiny elements on the diagonal of R? 
+    % If so, it means that some columns of A are linearly dependent.
+    zeroDiagIndex = findSmallIndices(diag(R));
+
+    % If there are fewer than DIM+1 ellipses, "add a fictional zero row to the end of R",´.
+    % This means that the (fictional) last element of the diagonal is zero 
+    % (i.e. it's index belongs in the list of small elements).
+    % This will trigger the correct solution option later. Don't worry, no panic, read on.
+    if isempty(zeroDiagIndex) &&  numFoci <= DIM
+        zeroDiagIndex = DIM + 1;
+    end
+    
+    % If there are more than one zero on the diagonal, the intersections can not be solved.
+    if length(zeroDiagIndex) > 1
+        return 
+    end
+
+    % If there are no zeros on the diagonal, the problem can be solved
+    if isempty(zeroDiagIndex)
+        xyw = T \ r; % == R \ (Q' * r);
+        intersectionPoints = xyw;
+        return
+    end
+
+    % Now matrix T is made of 3 columns: T = [A B C] and the problem to be solved is
+    % T * [x y w]' = Ax + By + Cw = r  such that  x^2 + y^2 = w^2
+    % Now one of the diagonals of the upper diagonal matrix R is zero,
+    % meaning that the corresponding column of T can be expressed as a linear combination
+    % of the preceeding columns. 
+    % Let's go over all possible positions of the zero diagonal and see what it means.
+    switch zeroDiagIndex
+        case 1
+            % R(1,1) = 0, meaning that column A = 0. 
+            % Solve By + Cw = r and use constraint x^2 = w^2 - y^2
+            yw = T(:,2:3)\r; % Now yw(1) = y, yw(2) = w
+            x = sqrt(yw(2)^2 - yw(1)^2);
+            x = [x  -x];
+            y = [yw(1) yw(1)];
+            w = [yw(2) yw(2)]; % Uncomment if w is needed
+        case 2
+            % R(2,2) = 0, meaning that column B = g * A for some g.
+            % Let's find what g is:
+            g = R(1,2) / R(1,1);
+            % Now we have  Ax + (g * A)y + Cw = r  such that  x^2 + y^2 = w^2
+            % --> A (x + g*y) + Cw = r. 
+            % Denote h=x + g*y, meaning x = h - g*y.
+            % We can now solve h and w without using column B:
+            hw = T(:,[1 3]) \ r; % Now hw(1) = h, hw(2) = w
+            % Now w^2 = x^2 + y^2 = (h - g*y)^2 + y^2
+            % --> (g^2 + 1)*y^2 - 2*g*h*y + (h^2 - w^2) = 0.
+            %  This is a 2nd degree polynomial on y. Let's solve it.
+            y = poly2zeros(g^2 + 1, -2*g*hw(1), hw(1)^2 - hw(2)^2);
+            % Now that we know y, we also know x = -g*y + h
+            x = hw(1) - g*y ;
+            w = [hw(2) hw(2)]; % Uncomment if w is needed
+        case 3
+            % R(3,3) = 0, meaning that column C = g1*A + g2*B for some g1,g2.
+            % Let's find what g is:
+            g = R(1:2,1:2) \ R(1:2,3);
+            % Now A*x + B*y + (g1*A + g2*B)*w = r such that  x^2 + y^2 = w^2
+            % --> A*(x+g1*w) + B*(y+g2*w) = r.
+            % Denote h1 = x+g1*w, h2 = y+g2*w, meaning x = h1 - g1*w, y = h2 - g2*w
+            % Now solve h1 and h2 using columns A and B
+            h = T(:,1:2) \ r;
+            % Now  w^2 = x^2 + y^2 = (h1 - g1*w)^2 + (h2 - g2*w)^2
+            % --> (g1^2 + g2^2 - 1)*w^2 - 2*(g1*h1+g2*h2)*w + (h1^2+h2^2) = 0
+            %  This is a 2nd degree polynomial on w. Let's solve it.
+            w = poly2zeros(g(1)^2 + g(2)^2 - 1, -2*(g(1)*h(1) + g(2)*h(2)), h(1)^2 + h(2)^2);
+            % Now we know w, g and h so we know x and y
+            xy = h - g*w;
+            x = xy(1,:);
+            y = xy(2,:);
+        otherwise
+            error("small diags should be 1 2 or 3 (" + zeroDiagIndex + ")")
+    end
+    intersectionPoints = [x;y;w];
+end
+
+function intersectionPoints = solveSpheroidIntersections3D(fp, diam)
+%   Inputs:
+%       fp - Second focal points as a 3-by-N matrix (the first focal points are always at the origin) 
+%       diam - Diameters of the ellipses along the major semiaxes, N-by-1 vector.
+%   Output:
+%       intersectionPoints - Coordinates of the points where all N ellipses intersect, or an empty vector if they don't intersect.
+    DIM = 3;
+    intersectionPoints = [];
+    numFoci = length(diam); % Number of ellipses (must be DIM or DIM+1)
+
+%  An intersection point (x,y,z) of the spheroids and it's distance w = sqrt(x^2 + y^2 + z^2) from the origin 
+% (i.e. the second focal point) satifies the following conditions:
+%       (x - fp(1,1))^2 + (y - fp(2,1))^2 + (z - fp(3,1))^2 = (diam(1) - w)^2
+%       ...
+%       (x - fp(1,N))^2 + (y - fp(2,N))^2 + (z - fp(3,N))^2 = (diam(N) - w)^2
+%                                           x^2 + y^2 + z^2 = w^2
+% By expanding the squares in N first equations and using the last equation, 
+% we get a linear problem T * [x y z w]' = r  with matrix T and vector r defined as as follows:
+    T = -2 * [fp(1,:)' fp(2,:)' fp(3,:)' diam];
+    r = diam.^2 - fp(1,:)'.^2 - fp(2,:)'.^2 - fp(3,:)'.^2;
+    
+    % QR transform will later help with sorting out corner cases later
+    [~,R] = qr(T);
+    
+    % Are there tiny elements on the diagonal of R? 
+    % If so, it means that some columns of A are linearly dependent.
+    zeroDiagIndex = findSmallIndices(diag(R));
+
+    % If there are fewer than DIM+1 ellipses, "add a fictional zero row to the end of R",´.
+    % This means that the (fictional) last element of the diagonal is zero 
+    % (i.e. it's index belongs in the list of small elements).
+    % This will trigger the correct solution option later. Don't worry, no panic, read on.
+    if isempty(zeroDiagIndex) && numFoci <= DIM
+        zeroDiagIndex = DIM + 1;
+    end
+    
+    % If there are more than one zero on the diagonal, the intersections can not be solved.
+    if length(zeroDiagIndex) > 1
+        return 
+    end
+    
+    % If there are no zeros on the diagonal, the problem can be solved
+    if isempty(zeroDiagIndex)
+        xyzw = T \ r; % == R \ (Q' * r);
+        intersectionPoints = xyzw;
+        return
+    end
+    % Now matrix T is made of 4 columns: T = [A B C D] and the problem to be solved is
+    % T * [x y z w]' = Ax + By + Cz + Dw = r  such that  x^2 + y^2 + z^2 = w^2.
+    % Now one of the diagonals of the upper diagonal matrix R is zero,
+    % meaning that the corresponding column of T can be expressed as a linear combination
+    % of the preceeding columns. 
+    % Let's go over all possible positions of the zero diagonal and see what it means.
+    switch zeroDiagIndex
+        case 1
+            % R(1,1) = 0, meaning that column A = 0. 
+            % Solve By + Cz + Dw = r and use constraint x^2 = w^2 - y^2 - z^2
+            yzw = T(:,2:4)\r; % Now yzw(1) = y, yzw(2) = z, yzw(3) = w
+            x = sqrt(yzw(3)^2 - yzw(1)^2 - yzw(2)^2);
+            x = [x  -x]; % x and -x both are solutions.
+            y = [yzw(1) yzw(1)];
+            z = [yzw(2) yzw(2)];
+            w = [yzw(3) yzw(3)]; % Uncomment if w is needed
+        case 2
+            % R(2,2) = 0, meaning that column B = g * A for some g.
+            % Let's find what g is:
+            g = R(1,2) / R(1,1);
+            % Now we have  Ax + (g * A)y + Cz + Dw = r  such that  x^2 + y^2 + z^2 = w^2
+            % --> A (x + g*y) + Cz + Dw = r. 
+            % Denote h=x + g*y, meaning x = -g*y + h.
+            % We can now solve h and w without using column B:
+            hzw = T(:,[1 3 4]) \ r; % Now hzw(1) = h, hzw(2) = z, hzw(3) = w
+            % Now w^2 = x^2 + y^2 + z^2 = (-g*y + h)^2 + y^2 + z^2
+            % --> (g^2 + 1)*y^2 - 2*g*h*y + (h^2 + z^2 - w^2) = 0.
+            %  This is a 2nd degree polynomial in y. Let's solve it.
+            y = poly2zeros(g^2 + 1, -2*g*hzw(1), hzw(1)^2 + hzw(2)^2 - hzw(3)^2);
+            % Now that we know y, we also know x = -g*y + h
+            x = -g*y + hzw(1);
+            z = [hzw(2) hzw(2)];
+            w = [hzw(3) hzw(3)]; % Uncomment if w is needed
+        case 3
+            % R(3,3) = 0, meaning that column C = g1*A + g2*B for some g1,g2.
+            % Let's find what g is:
+            g = R(1:2,1:2) \ R(1:2,3);
+            % Now A*x + B*y + (g1*A + g2*B)*z + Dw = r such that  x^2 + y^2 + z^2 = w^2
+            % --> A*(x+g1*z) + B*(y+g2*z) + Dw = r.
+            % Denote h1 = x+g1*z, h2 = y+g2*z, meaning x = h1 - g1*z, y = h2 - g2*z
+            % Now solve h1 and h2 using columns A, B and D
+            hw = T(:,[1 2 4]) \ r; % Now h(1:2) are h1,h2 and hw(3) = w
+            % Now w^2 = x^2 + y^2 + z^2 = (h1 - g1*z)^2 + (h2 - g2*z)^2 + z^2
+            % --> (g1^2 + g2^2 + 1)*z^2 - 2*(g1*h1 + g2*h2)*z + (h1^2 + h2^2 - w^2) = 0
+            % This is a 2nd degree polynomial in z. Let's solve it
+            z = poly2zeros(g(1)^2 + g(2)^2 + 1, -2*(g(1)*hw(1) + g(2)*hw(2)), hw(1)^2 + hw(2)^2 - hw(3)^2);
+            xy = hw(1:2) - g*z;
+            x = [xy(1,1) xy(1,2)];
+            y = [xy(2,1) xy(2,2)];
+            w = [hw(3) hw(3)];
+        case 4
+            % R(4,4) = 0, meaning that column D = g1*A + g2*B + g3*C for some g1,g2,g3.
+            % Let's find what g is:
+            g = R(1:3,1:3) \ R(1:3,4);
+            % Now A*x + B*y + C*x + (g1*A + g2*B + g3*C)*w = r such that  x^2 + y^2 +z^2 = w^2
+            % --> A*(x+g1*w) + B*(y+g2*w) + C*(z+g3*w) = r.
+            % Denote h1 = x+g1*w, h2 = y+g2*w, h3 = z+g3*w, 
+            % meaning x = h1 - g1*w, y = h2 - g2*w, z = h3 - g3*w
+            % Now solve h1, h2, h3 using columns A, B and C
+            h = T(:,1:3) \ r;
+
+            % Now w^2 = x^2 + y^2 + z^2 = (h1 - g1*w)^2 + (h2 - g2*w)^2 + (h3 - g3*w)^2
+            % --> (g1^2 + g2^2 + g3^2 - 1)*w^2 - 2*(g1*h1+g2*h2+g3*h3)*w + (h1^2+h2^2+h3^2) = 0
+            %  This is a 2nd degree polynomial on w. Let's solve it.
+            w = poly2zeros(g(1)^2 + g(2)^2 + g(3)^2 - 1, ...
+                           -2*(g(1)*h(1)+g(2)*h(2)+g(3)*h(3)), ...
+                           h(1)^2 + h(2)^2 + h(3)^2);
+            % Now we know w, g and h so we know x, y and z
+            xyz = h - g*w;
+            x = xyz(1,:);
+            y = xyz(2,:);
+            z = xyz(3,:);
+        otherwise
+            error("small diags should be 1 2 or 3 (" + zeroDiagIndex + ")")
+    end
+    intersectionPoints = [x;y;z;w];
+end
+
+function z = poly2zeros(a,b,c)
+% poly2zeros  Finds zeros of polynomial a*x^2 + b*x + c
+    determinant = sqrt(b*b - 4*a*c);
+    z = (-b + [determinant -determinant]) / (2*a);
+end
+
+function smallIndices = findSmallIndices(v)
+% smallIndidices Finds "very small" indices of vector v
+  smallIndices = find(abs(v) <= sqrt(eps(norm(v))));
+end
diff --git a/octave/threePointDistance.m b/octave/threePointDistance.m
new file mode 100644
index 0000000..02cd776
--- /dev/null
+++ b/octave/threePointDistance.m
@@ -0,0 +1,11 @@
+function distance = threePointDistance(point1Array, point2, point3)
+%measuredDistances Calculate distance from each Point1 to Point2 and further on to Point3 
+% Inputs:
+%   point1Array - array of Point1 coordinates, DIM x N array
+%   point2 -  DIM x 1 array
+%   point3 - DIM x 1 array
+% Output:
+%  distances - point1(k) -> point2 -> point3 distance, 1 x N array
+    distance = vecnorm(point1Array - point2, 2, 1);
+    distance = distance + norm(point2 - point3);
+end
\ No newline at end of file
diff --git a/pictures/2-ellipses-case-1.png b/pictures/2-ellipses-case-1.png
new file mode 100644
index 0000000..31562ba
Binary files /dev/null and b/pictures/2-ellipses-case-1.png differ
diff --git a/pictures/2-ellipses-case-2.png b/pictures/2-ellipses-case-2.png
new file mode 100644
index 0000000..d8d4e14
Binary files /dev/null and b/pictures/2-ellipses-case-2.png differ
diff --git a/pictures/2-ellipses-case-3.png b/pictures/2-ellipses-case-3.png
new file mode 100644
index 0000000..6a28abf
Binary files /dev/null and b/pictures/2-ellipses-case-3.png differ
diff --git a/pictures/3-ellipses-case-1.png b/pictures/3-ellipses-case-1.png
new file mode 100644
index 0000000..dc43e39
Binary files /dev/null and b/pictures/3-ellipses-case-1.png differ
diff --git a/pictures/3-ellipses-case-2.png b/pictures/3-ellipses-case-2.png
new file mode 100644
index 0000000..a1305c6
Binary files /dev/null and b/pictures/3-ellipses-case-2.png differ
diff --git a/pictures/3-ellipses-case-3.png b/pictures/3-ellipses-case-3.png
new file mode 100644
index 0000000..f8892ab
Binary files /dev/null and b/pictures/3-ellipses-case-3.png differ
diff --git a/pictures/3-ellipses-case-4.png b/pictures/3-ellipses-case-4.png
new file mode 100644
index 0000000..780b6fc
Binary files /dev/null and b/pictures/3-ellipses-case-4.png differ
diff --git a/pictures/3-spheroids-case-1.png b/pictures/3-spheroids-case-1.png
new file mode 100644
index 0000000..14dfe32
Binary files /dev/null and b/pictures/3-spheroids-case-1.png differ
diff --git a/pictures/3-spheroids-case-2.png b/pictures/3-spheroids-case-2.png
new file mode 100644
index 0000000..1b4499b
Binary files /dev/null and b/pictures/3-spheroids-case-2.png differ
diff --git a/pictures/3-spheroids-case-3.png b/pictures/3-spheroids-case-3.png
new file mode 100644
index 0000000..9027849
Binary files /dev/null and b/pictures/3-spheroids-case-3.png differ
diff --git a/pictures/3-spheroids-case-4.png b/pictures/3-spheroids-case-4.png
new file mode 100644
index 0000000..836f758
Binary files /dev/null and b/pictures/3-spheroids-case-4.png differ
diff --git a/pictures/4-spheroids-case-1.png b/pictures/4-spheroids-case-1.png
new file mode 100644
index 0000000..4767cff
Binary files /dev/null and b/pictures/4-spheroids-case-1.png differ
diff --git a/pictures/4-spheroids-case-2.png b/pictures/4-spheroids-case-2.png
new file mode 100644
index 0000000..941cfcb
Binary files /dev/null and b/pictures/4-spheroids-case-2.png differ
diff --git a/pictures/4-spheroids-case-3.png b/pictures/4-spheroids-case-3.png
new file mode 100644
index 0000000..b5b0131
Binary files /dev/null and b/pictures/4-spheroids-case-3.png differ
diff --git a/pictures/4-spheroids-case-4.png b/pictures/4-spheroids-case-4.png
new file mode 100644
index 0000000..10ccc7c
Binary files /dev/null and b/pictures/4-spheroids-case-4.png differ
diff --git a/pictures/ellipse-demo.png b/pictures/ellipse-demo.png
new file mode 100644
index 0000000..d238b45
Binary files /dev/null and b/pictures/ellipse-demo.png differ
diff --git a/pictures/example-ellipse-2.png b/pictures/example-ellipse-2.png
new file mode 100644
index 0000000..30bfb77
Binary files /dev/null and b/pictures/example-ellipse-2.png differ
diff --git a/pictures/example-ellipse-3.png b/pictures/example-ellipse-3.png
new file mode 100644
index 0000000..8cd2a8f
Binary files /dev/null and b/pictures/example-ellipse-3.png differ