As machine learning algorithms become increasingly complex, the need for transparent and interpretable models grows more critical. Shapley values, a local explanation method derived from cooperative game theory, is an explanatory method that describes the feature attribution of machine learning models. This is defined by the contribution a feature has to one single prediction. In this thesis the Shapley value formula is detailed, along with its properties. Three different ways to approximate the Shapley values are then introduced: Monte Carlo method for value function, Monte Carlo method through permutations and Kernel SHAP. Implementation of these methods using different datasets reveals that while both methods produce nearly identical Shapley values, Kernel SHAP is significantly faster. This research contributes to the field by demonstrating the advantages of integrating Shapley values into machine learning models, to enhance model transparency and trustworthiness.

Gross Primary Production (GPP) is the amount of carbon dioxide (CO2) that is fixed by an ecosystem through photosynthesis. It is a key variable to understand the carbon cycle and biodiversity conservation. In-situ measurements of GPP are possible at a local scale. To estimate GPP on a global scale, models based on remote sensing data are used. One of them is the Penman-Montheith-Leuning (PML) model, which is constructed on physical processes.

In this thesis, a sensitivity analysis of the PML model is performed for the case study of Torgnon subalpine grassland. It is defined around three aspects: the application of the PML model to Torgnon subalpine grassland, the impact of the input data on the output data, and the reliability of the calibration procedure. We found out that the PML model can be applied to the study site, but the value of the parameters presented in the literature are not optimal (R2 = 0.57 using the literature calibration and R2 = 0.94 in the optimised case).

Moreover, the model is extremely sensitive to the Leaf Area Index (LAI), and, to a lesser extent, to temperature and Photosynthetically Active Radiations (PAR). However, it is weakly dependent on the CO2 concentration of the atmosphere. The analysis of the calibration procedure showed that a simplified PML model, that does not consider the CO2 concentration and Vapour Pressure Deficit (VPD), performs similarly to the initial model. An ensemble estimation of GPP is obtained by using an ensemble of optimal parameters. It suggests that errors in the calibration principally impact the spring and summer estimations.

This paper primarily employs Item Response Theory (IRT) to estimate item characteristics and the proficiency levels of students as reflected in the exam results. The process includes the application of algorithms for item characteristic parameter estimation and the utilization of statistical techniques for item selection from an item pool. Additionally, Classical Test Theory (CTT) is also used to gain insights into the item characteristics. The mathematical frameworks behind every algorithm and model will be introduced in detail. Our ultimate objective is to create an item bank that unifies all items from various exam versions onto a common scale. The technique to put items from different versions on the same scale is the test equating technique, and it will also be described in detail.

After a comprehensive analysis to determine the appropriate model within the framework of Item Response Theory (IRT), the Rasch model has been selected for all versions of the exams. One item from version 9 has been removed from the item pool due to its unsatisfactory item fit. Based on the Item-Map plot, the conclusion can be made that all three exam versions appear to be relatively easy for students to answer, as evidenced by the high bar of the exams. Ultimately, the item bank has been successfully established through the application of two test-equating methods, with results indicating its reliability.

We introduce a method of designing NMR pulse sequences, with a specific focus toward complete decoupling of carbon-13 spin qubits coupled to a nitrogen-vacancy (NV) center in diamond. Using Average Hamiltonian Theory, we calculate different intermediate Hamiltonians corresponding to different rotations implemented by NMR pulses. Two novel pulse sequences are obtained containing 24 and 12 evolution periods respectively. Assuming ideal and instantaneous pulses, performance for both sequences is better than the WAHUHA + echo sequence for total evolution times of less than 5 ms and comparable for greater evolution times. Analysis of sensitivity points toward non-robustness for angle errors in the applied pulses when the direction of pulses is not taken into account, however this can be corrected for using chirality sums. Suggestions for further research include the study of non-globally applied pulses and application of the design method to non-zero effective Hamiltonians.

Low field magnetic resonance imaging (MRI) scanners provide a unique low-cost alternative to conventional MRI scanners. Nevertheless, low-field scanners come with drawbacks such as reduced signal-to-noise ratio and resolution, and also distorted images caused by magnetic field inhomogeneity and non-linear gradient fields. Despite this, it still provides a more accessible way to provide MRI in resource-limited areas. The main goal of this thesis is to develop an algorithm that can reconstruct 3D data from the low-field scanner efficiently and without distortion to the image. To do this, conjugate phase reconstruction (CPR) is employed, particularly frequency segmented reconstruction (Noll, 1991) and multifrequency interpolation (Man et al., 1997).

This thesis contains the development of continuous Kepler orbit- and a discrete numerical integration-based collision detection algorithms in a system of LEO satellites, which in combination with collision algorithm form a simplified space debris evolution model. This model is then used to study the Kessler syndrome. The continuous and discrete algorithms get their names from the solutions of the Two Body Problem (TBP) and the methods for collision detection that they are based on; the analytical and continuous time solution of TBP resulting in the Kepler orbits and the numerical, discrete time Velocity Verlet integration of the TBP. The collision model consists of an algorithm for fragmentation collisions largely based on the NASA Standard Breakup Model and a method for elastic, random scattering collisions. Comparison between the continuous and discrete algorithms shows that on average both predict the same time to the first collision in a system of homogeneously distributed satellites. The algorithms differ in their efficiency depending on the number and the radius of the satellites in and the geometry of the system. For relatively small satellite numbers in large systems, the continuous algorithm is computationally more efficient. However, as more satellites or fragments result from previous collision, the continuous algorithm is outperformed by the discrete algorithm. Consequentially, its time complexity appears to be O(N²). Armed with this knowledge, the continuous algorithm is used to show that an initially small system of satellites is able to evolve into a large population of debris particles within several decades. Similarly, the discrete algorithm is used to show that an ordered collection of satellites in an homogeneously distributed system of debris-like particles exhibits the effect that a collision early on in the simulation can cause a cascade of collisions at a later stage. Hence Both the discrete and continuous algorithms predict a Kessler Syndrome and mimic predictions made by more advanced models from leading space agencies like NASA’s LEGEND, ESA’s DELTA and JAXA’s LEODEEM [Lio+13].Future research could focus on including atmospheric drag and gravitational perturbations to the continuous algorithm, thereby lengthening the time frame during which it can realistically simulate a system of satellites in LEO. To achieve this, it is suggested that one execute the calculations inherent to the algorithm in parallel on a GPU, as these are independent of each other. 

Emission of ammonia and nitrogen oxides puts pressure on vulnerable ecosystems in The Netherlands. Satellite data from TROPOMI and CrIS gives vertical column density (VCD) maps for nitrogen dioxide and ammonia respectively. The Flux-Divergence method converts a VCD map to an emission map, by temporally averaging the divergence of the flux of the trace gas and estimating a sink term. This research proposes two improvements to the current implementation of the Flux-Divergence method at the KNMI for nitrogen dioxide. On the one hand, noise is reduced by computing the divergence before interpolating TROPOMI data to a regular grid. On the other hand, some emission locations are enhanced by reducing divergence in the wind field. An algorithm is provided to reduce this divergence. This thesis also explores the sensitivities of the method to different choices in its implementation. Finally, this research implements the Flux-Divergence method for ammonia. The thesis shows that, for ammonia, the method depends strongly on the estimation of the sink term, as the flux-divergence map is not able to capture homogeneous, spread-out emission sources.

Ammonia (NH3) is an important chemical compound in the nitrogen cycle. Ammonia is an essential nutrient and an important part of fertilizer, which in soil leads to increased growth of crops. However, the current excess of NH3 emissions is a hazard to environmental and human health. While ammonia
emissions need to be reduced, the emission estimates are still highly uncertain.
In this study, a method is developed to combine the chemical transport model LOTOS-EUROS with measured ammonia concentrations to improve the ammonia emission estimates in the Netherlands and the surrounding regions. The measurements used are generated by the miniDOAS instruments on LML stations and the IRS instrument on board the future MTG-S satellite. The proposed method is an adjoint-free 4DVar method. The 4DVar method aims to retrieve the emission parameters for which the LOTOS-EUROS model determines NH3 concentrations that resemble the measurements while keeping the emissions fairly similar to the original inventories. A linear approximate model has been developed, which uses the near-linearity of the NH3 concentrations in terms
of the NH3 emissions. When using the approximate model, the 4DVar
method can be solved without using an adjoint model, making the method adjoint-free. Subsequently, the 4DVar cost function is rewritten to allow the emission parameters to have a lognormal prior distribution. A maximum likelihood approach is developed to estimate the parameters of both the prior distribution of the emissions and the likelihood of the measured observations. Last, a preconditioner in reduced space has been considered, to estimate emissions on a fine spatial resolution, while keeping the computational cost
feasible. This preconditioner uses the property that the emission parameters are correlated in space. First, the methodology has been tested in an identical twin experiment where the emissions vary only in time and strength, using the LML observations. It was concluded that the methodology worked well
for short periods (less than 30 days), but the results were dominated by observational noise. When using the real observations in the 4DVar
method, the results seem unrealistic Second, the method has been tested in an identical twin experiment where emissions vary in space, as well as in time and strength. Here, synthetic IRS observations of the future MTG-S satellite are
used. The optimized emissions did resemble the true emissions of the twin experiment. Observational noise appeared to no longer be an issue. However, the results were not perfect. The regions with the highest emission increase appeared to be underestimated, low emission areas appeared to have large
relative estimation errors, and estimates for coastal regions seemed to be incorrect. Hence, on a local scale, the emission estimates can be imperfect, but overall the adjoint-free 4DVar method does greatly improve the emission inventories.

A severe incident with the cargo ship MSC Zoe, which lost hundreds of containers near the Dutch coast, shows how important it can be to detect and localise hidden objects at sea quickly. Also in defence applications localisation of hidden objects in water is important. If the objects are made of magnetic materials, they disturb the Earth magnetic field. This disturbance is called the object's magnetic signature. An aerial drone that carries magnetic field sensors could in principle be able to locate magnetic objects autonomously. Another defence-related application of such a drone would be estimating magnetic signatures of ships under water, using measurements above water - signature translation. For both applications, algorithms need to be developed that process magnetic field data, and computes estimates. In this thesis, data-driven techniques are applied to the localisation problem and the magnetic signature translation problem.

For the localisation problem, an algorithm is proposed to localise a magnetic dipole using a limited number of noisy measurements from a sensor array forming a horizontal grid. The algorithm is based on the theory of compressed sensing. In the algorithm, a number of sensors is chosen which perform each a measurement of the three magnetic field components. The sensors can be chosen randomly from the sensor array, but also optimal placement using QR pivoting is considered. Using the obtained field measurements, a sparse representation in the location domain is computed using L1-optimisation. Based on the resulting sparse representation, a classification is produced, which consists of estimates on its location and on its magnetic moment magnitude and orientation. The possibility of performing iterations is explored, where the basis and chosen sensors improve after every location estimate. Using results from simulations as well as experiments, the algorithm is shown to be effective in localising magnetic dipoles.

To solve the signature translation problem, two approaches are investigated. One algorithm was developed using the Gappy POD technique, and one using neural networks. Using only a few measurements (a gappy measurement), a full signal can be reconstructed. This method is adapted to also be able to translate from a field above a ship to a field below a ship. Tikhonov regularisation is applied in the process to obtain better results. Second, different neural networks are created, trained using data above and below a ship. Linear and non-linear networks are compared, and standard loss functions are compared to physics-guided losses. Both the Gappy POD based algorithm and linear neural networks are shown to give good magnetic signature estimates. The Gappy POD method is greatly improved by applying Tikhonov regularisation. Results show to improve when more basis modes are considered, and when more sensors are used for measurements. Linear neural networks show the same results as Gappy POD using Tikhonov regularisation, indicating that some regularisation is done while training the neural network. Non-linear neural networks show no improvement over linear ones, and a physics-guided loss function is not needed for a resulting field that obeys known laws of magnetism.

Magnetic resonance imaging (MRI) scanners are a crucial diagnostic tool for radiologists. They are able to render two­ and three­dimensional images of the body without exposure to harmful radiation. MRI systems are, however, costly to build and maintain. This adversely impacts access to these scanners in developing regions. In an effort to combat this problem, a low­field MRI scanner is being developed. Conventional MRI scanners utilize a superconducting solenoid to generate the main magnetic field. The low­field scanner, on the other hand, induces the main magnetic field through a Hallbach array of permanent neodymium magnets. While beneficial for production and maintenance costs, as well as portability, the Hallbach array is not able to generate a perfectly homogeneous magnetic field. The inhomogeneities present in the main magnetic field result in distortion of the images when reconstructed using conventional fast Fourier transform (FFT) methods. To counteract this, a reconstruction method that utilizes field information needs to be employed. In this thesis, existing methods to determine and utilize the field information to correct image distortion are explored. From this analysis, it is evident that model­based (MB) methods are most suitable for reconstruction of data from the low­field scanner. Current MB methods are only implemented for two­dimensional reconstruction. The goal of this thesis is to expand these methods to three­dimensional reconstruction. A novel MB method for three­dimensional reconstruction is presented. This new method is able to circumvent memory constraints that arise from reconstruction of large data sets. Though the new method requires several hours to reconstruct a 128 × 128 × 30 data set, visual inspection indicates that an accurate result is achieved.

In recent years many automatic differentiable programming frameworks have been developed in which numerical programs can be differentiated through automatic differentiation (AD). Examples of these frameworks are Theano, TensorFlow and Pytorch. These frameworks are widely used in Machine Learning. AD also finds applications in the field of computational fluid dynamics (CFD). It is used to develop discrete adjoint CFD code for research concerning for instance sensitivity analysis, data assimilation and design optimization. However, the use of the automatic differentiable programming frameworks in the field of CFD is limited. One can find some examples in the literature on how to find a numerical solution to an initial value problem using a differentiable programming framework. In this work it will be clarified how one can implement an semiimplicit
time integration scheme for a staggered grid to simulate the propagation of long waves in water with a free surface in TensorFlow. A main advantage of the automatic differentiable programming frameworks is the user friendly application
programming interface (API) for AD. No research has been conducted to use this API in the field of CFD. In this work an example will be given how one can use TensorFlow for research concerning sensitivity analysis. AD requires a significant allocation of memory on a CPU/GPU when working with fine meshes and/or long simulations and since CPU/GPU memory is finite, the method checkpointing is
proposed to make it feasible to perform sensitivity analysis when working with fine meshes and/or long simulations. Another main advantage of the differentiable programming framework TensorFlow is the use of compute unified device architecture (CUDA) of a NVIDIA GPU in order to perform computations in parallel, which results in a significant reduction in computation time. A Benchmark will be given that indicates the computational efficiency of TensorFlow compared to a loop over grid implementation in NumPy and a Fortran CPU scalar implementation.

In the fight against money laundering, demand for data-driven Anti-Money Laundering (AML) solutions is growing. Particularly anomaly detection algorithms have proven effective in the detection of suspicious customer behaviour, as well as observing patterns otherwise hidden in customer transaction data. In this thesis, the Isolation Forest anomaly detection algorithm is studied in combination with the model-specific local explanation method, Multiple Indicator Local Depth-based Isolation Forest Feature Importance (MI-Local-DIFFI). To expand Isolation Forest to mixed-attribute data sets, the incorporation of nominal features is explored in more detail. This analysis resulted in the introduction of Isolation Forest with Categorical Sampling (iForestCS ), a methodology that directly incorporates nominal attributes into an isolation tree without the need of encoding it onto a numerical scale. This method is tested against different encoding strategies and Isolation Forest Conditional Anomaly Detection (iForestCAD) using different synthetic data sets. The method shows improved performance to the utilization of encoding strategies for different parameters of the underlying synthetic data. Furthermore, this thesis explores the potential of ternary Isolation Forest, in which the branching strategy of an isolation tree is expanded to produce three child nodes. It is demonstrated using synthetic data, that particularly the performance of MI-Local-DIFFI reduces when applied to a ternary Isolation Forest. Finally, the research considers a practical use-case. Using customer transaction data from Triodos Bank, the locally explainable Isolation Forest is applied to mixed-attribute customer transaction data. This has provided useful insight and resulted in the detection of suspicious customer behaviour and the introduction of new rules into business practices. Although the most interesting customer behaviour did not directly emanate from the nominal attributes, the method of incorporating nominal features resulted in differences when considering the anomalies with the highest anomaly scores.

Magnetic Resonance Imaging is one of the most widely used imaging modalities nowadays and it performs especially well imaging human organs such as the brain and liver. One of its main limitations is the relatively long imaging times, to overcome this issue and speed up the data acquisition, several techniques such as Parallel Imaging or PI have been developed. These techniques require advanced hardware and software to be able to decrease the acquisition time. On the hardware side, a highly efficient insert gradient coil has been designed and built at the University Medical Center Utrecht. Specialized software has to be implemented to optimally make use of this hardware. One of the recently proposed PI methods called Wave-CAIPI has been proved to achieve a ninth fold acceleration factor without compromising image quality. This project aims to investigate the time gain that can be achieved when combing the insert gradient coil with a Wave-CAIPI strategy. Two main aspects are reviewed. The first one is the maximum achievable under-sampling factor that does not compromise image quality. The second one is the decrease in acquisition time that can be obtained when using the insert gradient coil compared to conventional gradient systems while maintaining image quality. To do so, the strategy has been implemented and extensive simulations have been performed to optimize the MR acquisition parameters. To prove the results from the simulations, the Wave-CAIPI sequence was implemented in a 7T scanner at the UMCU, where the acquired data was retrospectively under-sampled, obtaining the wave image to be further reconstructed. Limitations of previous works on Wave-CAIPI have been the gradient specifications, which can be overcome by the high-efficiency coil. It has been concluded that shorter acquisition times without compromising image quality are possible when using the insert coil compared with conventional systems. The time gain can be up to a factor of five, and sixteen fold under-sampling factors could be possible. The time gain can be especially useful for Echo Planar Imaging sequences, where switching faster gradients allows to acquire more signals in less time. The next steps for this research are to prospectively under-sample the data in a Wave-CAIPI fashion with a sixteen under-sampling factor and corroborate if sequences such as Echo Planar Imaging can be benefited with the time gain.

In dit project worden vier vereenvoudigingen van de trebuchet onderzocht. Er wordt gekeken naar het rendement met betrekking tot de energie overdracht van het zware tegengewicht op het afgeschoten voorwerp. De invloed van het toevoegen van extra armen en wielen wordt onderzocht. Bij meer energie overdracht kan het voorwerp verder weg geschoten worden en zal het voorwerp meer impact hebben als het op de grond terecht komt.

In financial and egineering problems, we are often faced with solving Partial-Integro Differential Equations (PIDEs). Rarely we can find an analytic solution in a closed form expression for these PIDEs, hence we turn to numerical schemes to accurately approximate the solution instead. Classically these methods are based on finite difference methods, however, we can turn certain kinds of PIDEs into a probabilistic representation, called Forward Backward Stochastic Differential Equations with Jumps (FBSDEJs). Solving the PIDE can now be done alternatively by solving a FBSDEJ. In this thesis we will first investigate the stochastic framework behind FBSDEJs and we will look into the uniqueness and existence of their solutions. Furthermore we propose a new numerical method which can efficiently solve FBSDEJs. The semi-discretisation is based on the classical Backward Differentiation Formula (BDF) methods, for the computation of the conditional expectations we use the COS method which makes use of Fourier cosine expansions, exploiting the knowledge we have about characteristic functions. Finally we implement the new method and we investigate it extensively both numerically and theoretically. We show that the BDFn schemes are highly stable and efficient for computing FBSDEJs, the initial steps still have to be investigated in greater detail so that we can make use of the high-order BDFn schemes.

In cases that the voltage source is made of multiple cos’s and we can compensate the current, it is always possible with the first Cauer form. This is caused by the shape of the found transfer function Z(s). This transfer function is a fraction where the odd powers of s are in the numerator, while in the denominator are only even powers of s. In this way the conditions 1 and 3 in chapter 4.2 are satisfied. For this special type of Z(s) satisfying condition 2, is having positive coefficients in the First Cauer form, hence the conclusion. In the tables in this thesis, it seems that when the first Cauer form works, the other forms areworking as well. It seems that the method used in this thesis to find a controller only works for small R, L and C. In the example the improved PF is optimal, like one would expect. I conclude that the method used in this thesis is working fine.

This thesis covers the power load flow problem and three numerical methods which can be used as a solution. We firstly discuss the concepts from electrical engineering required to discribe the problem. This includes alternating current and voltage and admittances of various electrical elements. We then use these concepts to derive a set of equations that govern how power flows through a network. To solve these equations, we consider the Newton-Raphson, line-search and trust-region algorithms. For these methods, we also look at its convergence and the amount of computational power required. We conclude in the cases where the parameters of the electrical network are within reasonable bounds, the Newton-Raphson algorithm is sufficient in precision and requires the least computational power. Otherwise, one of the other two methods may be tried.

The theory of morpho-poroelasticity is applied to tumor growth. This allows for the modeling of permanent deformation in the tissue as a result of the presence of tumor cells. The work done in this project can be divided into two parts. In the first part we review existing models for elasticity and build corresponding finite element models. The aim of the first part is to gain understanding in the behaviour of morphoelasticity and poroelasticity. In morphoelasticity the deformation tensor is decomposed into an elastic and plastic component, allowing for permanent deformations to occur in the tissue. The theory of morphoelasticity describes the interplay between a fluid in a porous tissue and the tissue's elastic properties. A common problem in poroelastic models is the occurence of non-physical oscillations in the pressure. The second part of this project contains novel work. First, a rigorous mathematical derivation is presented of the tuning parameter found in the diffusive stabilization method for poroelastic systems. Secondly, we present a finite element model for the novel combination of morpho- and poroelasticity. The derivation concerning the diffusive stabilization is also applicable to the morpho-poroelastic finite element model. We obtain a promising tumor growth model by combining a biochemcial tumor growth model with morpho-poroelasticity. This biochemical model is based on a nutrient transport equation and a tumor cell density evolution equation. Some example output of the tumor growth finite element model is shown. It can be used as a starting point for more elaborate tumor growth models.

In this report we present a model to simulate the performance and robustness of the Charge-Resonance Check (CR-Check) in the Nitrogen-Vacancy (NV) center. The CR-Check is a routine which verifies that the charge state of the NV center is NV$^-$ and that the lasers, addressing the NV center, are on resonance with the optical lines. This is done by shining a red (λ=637 nm) laser on the NV center and measuring the amount of fluorescence photons, since they are dependent on the resonance frequency being near the frequency of the laser. If these counts do not pass the threshold for success, then either a new measurement of the photons is done or the system is repumped with a green (λ=532 nm) high intensity laser if they are below the repump threshold. The most important parameters are the duration of the red laser (dt) and the spectral diffusion of the NV center (sigma), which is an intrinsic property of the NV center. The model was implemented in Python, where also Numba, a high-performance Python compiler was used. An algorithm for finding the optimal values is included, so that different implementations can be compared. Then parameter sweeps were performed to get insight on the effects of the parameters. The distribution of frequency spectrum was investigated, it follows the shape of the Generalized Normal distribution with a disturbance due to the memory parameter. We observed that the mean time till success goes up linearly with the sigma after it is higher than 80 MHz for dt. This effect can be explained by taking into account the Normal distribution which was used to determine the resonance frequency. It became clear that the N_thr_repump and N_thr_success often come as pairs, which have been shown to be linearly dependent on the rate. Also it has been observed that after a certain sigma amount the N_thr_repump and N_thr_success do not change anymore. We noticed that for fixed sigma there is a optimal value of dt. This optimal value of dt is dependent on sigma. Average number of photon counts was nearly independent of sigma and only linearly dependent on dt. This is a sign that the CR-Check can always bring the NV center sufficiently on resonance. It has been observed that when the NV center's spectral diffusion is below 40 MHz, it is significantly faster to use the CR-strategy with dt = 25 μs compared to the standard dt of 50 μs. Lastly new implementations, where p-values are used rather than the measured counts, were shown to be capable of performing a CR-Check procedure. These implementations have used either the cumulative density function (cdf) of the Poisson distribution or the t-test or the Wilcoxon signed-rank test. The test using the cdf performed very similar to the current implementation. The other two test were slower with regards to the mean time till success, as they seem to require a lower spectral diffusion than was mandated.

In this thesis the stability type of y=0 is being considered for the delay differential equation y''(t) + ay(t) + by(t-1) = 0 with a and b real numbers. It is already known that y=0 is stable when b=0 and a>0 and unstable when b=0 and a<0. The aim of this project is to determine the stability of y=0 for all values of a and b. First, the general stability theory for delay differential equations was highlighted before giving an in-depth stability analysis of the equation y''(t) + ay(t) + by(t-1) = 0. It turns out that a theorem of Pontryagin (1908 -1988) is really helpful for answering these stability questions. Due to this theorem all values for a and b are determined such that y=0 is asymptotically stable for y''(t) + ay(t) + by(t-1) = 0. However, this does not cover the stability type of y=0 for all values of a and b. So more analysis was done in order to give a full answer of the stability problem. The full answer was not achieved as there are still values for a and b where the stability is unknown. Finally, numerical solutions of y''(t) + ay(t) + by(t-1) = 0 are shown to confirm the results that are obtained.