Partial differential equations (PDEs) are of significant interest to scientists and engineers as they can be used to study a variety of physical and engineering problems such as transport processes (Bird et al., 2002) , structural mechanics (Zienkiewicz and Taylor, 2005) and electromagnetism, quantum physics (Gowers et al., 2008) to list a few. For engineers, PDEs are often one of the most complex mathematical topics that would be covered in the undergraduate course. However, solving PDEs for systems of academic and industrial interest often require numerical methods as analytical solutions may not be available. Therefore, subsequent courses on PDEs for engineers often focus on applications in different settings and employ numerical methods and/or empirical correlations to obtain usable results.
The importance of numerical methods and the development of solver codes is evident in the sheer diversity and sophistication of PDE solver packages and codes currently available. In addition, there is a large variety of ancilliary elements and tools such as pre and postprocessing softwares, wrappers and webapps with varying levels of functionalities. These codes and software are either available commercially or open-source or they can be developed in-house to meet specific needs. Educators face a complex challenge in deciding which set of tools should be used in their classroom. Multiple factors such as scope and objectives of the course, students' background and 1. Availability: Open-source software is much more available to users compared to commercial software which can be quite costly. 2. Ease of use: Factors such as the learning curve associated with using the software and/or the availability of reference material and support should be considered. 3. Scalability: Can the software be used for both smaller tasks such as testing and teaching and also be used for large problems in industrial and research settings. 4. Learning objectives: Is the class more theoretical and focused on the mathematics and physics of the problem and/or numerical method or is the focus on the application of a tool to solve an engineering problem. Another factor to consider is the academic background of the students e.g. is the course aimed at an undergraduate or graduate level.
To date, no code or software has been demonstrated to effectively address all four factors listed above in a teaching setting. Commercial software such as ANSYS and COMSOL are commonly used by industry and academics and scale very well. Furthermore, there is comprehensive reference material available and a variety of courses at the undergraduate and graduate level utilize these software, often in an application setting (Morscheidt et al., 2013; Vicéns and Zamora, 2014; Adair et al., 2014; Huang et al., 2019; Wicaksono et al., 2020) . However, these commercial packages fall short in terms of availability due to their high costs (Chen et al., 2014) and as many of the numerical complexities are abstracted away, they are not readily usable for teaching the fundamentals aspects of solving PDEs numerically.
MATLAB, whilst commercial in nature and therefore restricted in access, is often used by educators in numerical methods courses (Kiusalaas, 2015) . However, such student generated codes scale poorly as they are meant to be instructive in nature. These courses are also often at a senior undergraduate or graduate level which restricts the target audience. Educators also commonly use MATLAB and other similar commercially available software in application settings for problems involving systems of ordinary differential equations (ODEs) (Ibrahim, 2011; Li and Huang, 2017; Molina et al., 2019) or PDEs (Stephens and Mauzeroll, 2019) .
Fortunately, for many systems involving ODEs, including chemical kinetics, process modelling or complex systems, there are a wide variety of open-source libraries in Python such as SciPy (Virtanen et al., 2020) or other languages such as Julia (Rackauckas and Nie, 2017) which can be used effectively. These libraries and Python in general have found significant use in both research and teaching settings (Golman, 2016; Wang et al., 2017; Sawaki et al., 2020) .
However, the open-source landscape for solving PDEs is a bit more nuanced. Many open-source codes which are incredibly powerful and have been gaining traction in academia and industry such as OpenFOAM (Jasak et al., 2007) or FEniCS (Logg et al., 2012) have effectively addressed issues such as availability and scalability. However, care must be taken with regards to using these codes in the classroom as factors such as the ease of use may adversely impact meeting module learning objectives. For example, OpenFOAM is well known for having a steep learning curve (Chen et al., 2014) and is therefore illsuited for introductory courses on numerical methods or in an application setting.
We believe that the recent development of a variety of opensource PDE solver codes in Python presents an incredible opportunity for educators. The following list of solver codes is non-exhaustive and is categorised by discretisation method:
Finite volume method FiPy (Guyer et al., 2009) Finite element method FEniCS (Logg et al., 2012) , SfePy (Cimrman et al., 2019) Finite difference method Scikit FiniteDiff (Cellier and Ruyer-Quil, 2019) , Pystencils (Bauer et al., 2015) Spectral methods Dedalus (Burns et al., 2019) These solver codes can effectively address all four considerations. Being open-source means that these solvers are easily accessible to users. With the front-end written in Python, users are able to quickly understand and write code to solve problems of interest (Guyer et al., 2009 ). In addition, most of these solvers have active user discussion forums or means to get support from the developers and have comprehensive documentation available. Many of these codes are also actively being used in research and have been demonstrated to effectively scale (Logg et al., 2012; Burns et al., 2019; Cimrman et al., 2019) . Lastly, when collectively considered as a set of tools, these solver codes contain a wide variety of features and capabilities that can be used by users of varying skill level, giving educators a remarkable design space when considering integrating these codes into their course.
As a first iteration, we demonstrate the use of one of these codes, FiPy, as part of an introductory course to numerically solving PDEs with Python. A short course which ran over two days was developed and rolled out over the summer of 2020. We hope to integrate parts of this course into an advanced engineering mathematics module to complement the teaching of the PDE section of the course. The extensive course notes and all code used in the course will be made publicly available. We intend to supplement the code repository in due course with additional contributions to expand the material offered.
The course outlined in the present work is rolled out as a course run over two three-hour sessions held on two days in the summer of 2020. We intend to integrate components of the course into the PDE section of the second year Math II (Engineering Mathematics) module that students are required to take. The Math II module is quite comprehensive, covering topics such as vector calculus, Fourier series, PDEs and statistics. Prior to this module, students would have taken a first year mathematics module covering univariable calculus, linear algebra and ordinary differential equations along with an introductory coding course in MATLAB. Students would have also completed an introductory class on transport processes. Transport phenomena is chosen as the context for most of the exercises as students would already have a conceptual understanding of the problems. Educators in different disciplines can readily implement problems of interest to them with comparative ease.
To help students appreciate the whole pipeline of setting up a simulation and analysing the results, the course starts with a quick guide on how to set up Python and install FiPy. Four exercises of increasing complexity are developed for students to explore. Each exercise may be further broken down into parts to emphasise specific points of interest. Each exercise was designed to simultaneously introduce the necessary content at that stage and also lay the groundwork for students to progress to the next exercise to enable students to progress to more complex problems (Xie et al., 2020) . In exercise 1, students explore the Laplace equation in a problem from their module while in exercise 4, students consider the Cahn-Hilliard equation to model polymer blend demixing. Exercise 4 also introduces students to using ParaView (Ayachit, 2020) , a research-grade software often used by researchers to analyse and post-process simulation data.
As students may not have experience with Python and may be intimidated by the idea of numerically solving PDEs, much care was taken to ensure that the course is as accessible as possible. A detailed step-by-step guide is provided to set up the environment and install relevant software before the beginning of the course in a pre-course document. All code is written with extensive comments explaining what each line of code means. It is desired that students would see the code provided to them as a reference for them to use in the future. Furthermore, comprehensive notes with annotated screenshots were made for all parts of the course. The notes and other course material were made in collaboration with students to improve accessibility (Inguva et al., 2018; Xie et al., 2020) .
A pre and post course survey was designed and given to participants in the course to fill out anonymously. The survey was designed to evaluate students' perception of the course structure and approach and their self-efficacy and attitudes before and after the course. A series of statements were provided and students were asked to rate how strongly they agreed with the statements. A focus group discussion was also conducted to get students' feedback on how they perceived the course and also to garner feedback and suggestions on future strategies for course improvement and integration into the classroom.
We used the Spyder Integrated Developer Environment (IDE) which is distributed with the Windows Installation of Anaconda (Anaconda, 2016). Students are also given instructions on how to create an Anaconda environment as a tool for creating reproducible environments for simulations. The associated environment files for Windows and Linux installations are available in the code repository.
The first exercise serves as an introduction to working with FiPy and basic concepts in solving PDEs. The exercise aims to achieve the following:
1. Set up a simple case 2. Run the simulation 3. Visualise the results 4. Benchmark the results of the simulation to the analytical solution
The classic 2D Laplace equation on a unit square domain is used as the equation to be solved with the following boundary conditions:
where φ is the solution variable and x, y are the x and y coordinates, respectively. The analytical solution of the problem is written as follows:
φ(x, y) = ∞ n=0 4 sin((2n + 1)πx) sinh((2n + 1)πy) (2n + 1)π sinh((2n + 1)π) .
(3)
As shown in fig. 1 , the numerical and analytical solutions to the problem posed in exercise 1 are highly comparable. However, the impact of the truncation of the analytical solution can be seen whereby the solution at the top vertex is less smooth compared to the numerical solution. These figures are generated using the default FiPy viewer and students would be able to see these results as they execute the scripts in the IDE.
To help explain the results observed, the idea of steady-state conduction/diffusion on a plate was used to provide a physical intuition of the results. The vibrant colours of the result also helped facilitate introducing the idea of contour / isolines as lines of constant φ can be visually identified quite easily.
The second exercise introduces transient problems in the form of the familiar diffusion equation without or with a reaction. The exercise is broken down into three parts: In part 2A, the diffusion equation is solved with Dirichlet boundary conditions while in part 2B, Neumann boundary conditions are applied. In the last part 2C, a reaction term is added to the diffusion equation. All the problems are solved in 1D. Exercise 2 aims to achieve the following: The problem statements for each part are as follows:
where c is the concentration, t is the time and D is the diffusion coefficient which is set to 1.0. This problem has an analytical solution which can be written as follows:
The mesh Fourier number Fo Mesh is introduced as a stability criteria for explicit forward Euler time-stepping. As a starting point, students are provided a value of ∆t = 0.00018 which is stable.
Part 2B. The same equation eq. (4) as part 2A is solved, except the initial and boundary conditions are as follows: Part 2C. A reaction term is introduced to the diffusion equation as follows:
where κ is the rate constant and is set to −2.0. The initial and boundary conditions are the same as part 2A. The reaction term is treated as an implicit source term.
To fully appreciate the results from part 2A, students will need to explore a few settings namely the time-stepping scheme,∆t and the time-stride. A time-stride parameter is introduced into the code to allow students to visualise every nth step which is especially important when explicit time-stepping is used due to the small ∆t. When an explicit forward Euler time-stepping scheme is used, the default value of ∆t provided to students yields a stable solution. However, when ∆t = 0.0001998, numerical instabilities in the form of "wiggles" become noticeable as shown in fig. 2 and further increasing ∆t causes the unacceptably oscillatory behaviour which is the expected behaviour. However, using an implicit time-stepping scheme enables users to use much higher values of ∆t and still achieve stable results.
If students use a small time-step and a small time-stride (such as the given value of the time-step and a time-stride of 1), they would also be able to view the initial profile of the truncated analytical solution which has the effects of the higher frequency Fourier modes visible. As time progresses, the higher frequency modes decay very quickly and the truncated analytical solution corresponds to the numerical solution very well. This can help students to visualise various theoretical concepts pertaining to Fourier series.
Conceptually, parts 2B and 2C are quite similar. We present the expected results from part 2C in fig. 3 to demonstrate a typical plot students can generate if they use the TSVViewer in FiPy to export their data and generate more complicated plots on their own.
The CDR equation is of significant use to model the transport of chemical species. Part 3A of the exercise assumes no reaction and serves to visually demonstrate how the convective and diffusive transport function. Concepts such as the discretisation of the convection term and the derivation of the Courant number are covered in the notes. Part 3B considers a coupled system of equations with a non-linear reaction term. The aims of exercise 3 are as follows:
1. Visualise how convective and diffusive transport occur. 2. Appreciate the significance of Peclet number in convectiondiffusion problems 3. Appreciate the added complexity involved in modelling convection. 4. Appreciate how scaling a problem facilitates its solution 5. Implement a coupled equation system with a non-linear source term.
The problem statement for each part is as follows:
Part 3A. The convection-diffusion equation is first scaled and solved in non-dimensional form as follows: wheret is the dimensionless time and is scaled by the characteristic length L and velocity U as such t = L Ut .x is the dimensionless spatial coordinate and is scaled as such x = Lx. Lastly, Pe is the Peclet number and is given by Pe = UL D . U and L are set to 1.0 and Pe is set to 100. A 1D mesh of length 4 was generated with a mesh size ∆x = 0.01. The simulation is run until t = 2.0. Periodic boundary conditions are applied:
The concentration c is initialised with a narrow Gaussian distribution:
Part 3B. The simple chemical reaction between species A and B to form C which is shown as follows:
is given as first order with respect to both species A and B with the following simplified expression for the reaction rate:
where r i , i ∈ A, B, C is the reaction rate for each of the species A, B and C, κ is the rate constant and a, b, c is the concentration of species A, B and C, respectively. The equation system for the three species can be written as follows:
where U and D are the convective velocity and diffusion coefficient respectively and both are assumed to be constant and identical for all species. D, κ and U are set at 1.0, 30.0 and 3.0, respectively. The simulation domain is set at L = 4. Dirichlet boundary conditions are specified at the inlet as follows:
c(x = 0, t) = 0.0.
However, specifying the boundary conditions at the outlet is non-trivial. The appropriate boundary condition is termed the "free boundary condition", but is outside the scope of this course. Often, a zero-diffusive flux boundary condition is applied in such cases (Lynch et al., 2020) and it suffices for pedagogical purposes.
∂c ∂x (x = L, t) = 0.
This equation system introduces students to two interesting concepts in numerical methods. The first is the idea of strong vs weak coupling. The transport equations for species A and B are strongly coupled and need to be solved at the same time, however the transport equation for species C is weakly coupled and can be lagged in the solution process. The second concept is the idea of "sweeps". As the reaction terms are non-linear, the equation system is non-linear necessitating solving the linear equation system iteratively at that time-step using the result from the previous solution (Guyer et al., 2013) . Students may typically solve the steady-state version of this problem as a set of coupled ODEs numerically, but they now have the opportunity to explore the transient dynamics of the solution.
The effects of the different transport mechanisms are shown in fig. 4 . Students are able to solve each equation (convectiononly, diffusion-only and convection + diffusion) by commenting out a few lines in the code. The notes introduce the concept of numerical diffusion which is an important challenge in modelling convection and is demonstrated by the false diffusion present in the convection-only system as shown in fig. 4 . FiPy has multiple convection schemes present, and the VanLeer convection scheme was found to have the least amount of numerical diffusion and is therefore used. Students can explore how different parameters such as the mesh resolution and convection scheme can impact the results with minimal modification to the code. 
To wrap up the course, we wanted to expose students to a problem that is currently of significant academic interest and is comparatively simple to implement in FiPy. A modified Cahn-Hilliard equation system based on the work of Petrishcheva and Abart (2012) is selected as the model problem. This equation system can be used to study polymer blend demixing Nauman and He, 1994) and can be coupled with other equations to model processes such as membrane formation (Zhou and Powell, 2006) . The numerical solution is also challenging as the simulations tend to diverge at certain parameter values Brunswick et al., 1998) . The problem statement with derivation and assumptions is detailed in the course notes and students are given a brief outline of the problem itself with some examples of polymer blend demixing. While the following derivation is important to develop an understanding of how to simplify complex set of equations, students are not expected to follow this step by step. Emphasis is placed on a conceptual understanding of the Cahn-Hilliard equation and implementation and discussion of solution from FiPy as well as processing using ParaView.
We consider the case of a binary polymer blend consisting of species A and B. Polymers A and B have a chain length of N A and N B , respectively, and the interaction between the two species is characterised by the Flory-Huggins parameter χ AB . As a first set of parameters, students are provided the following values: N A = N B = 1000 and χ AB = 0.006. The initial volume fraction of species A, a 0 was set to 0.77 with random noise of ±0.03. The Cahn-Hilliard equation is a fourth order PDE that is typically split and solved as a set of coupled second order PDEs as such (Jokisaari et al., 2017) :
whereμ AB is the non-dimensional difference in chemical potentials between species A and B,t is the non-dimensional time, a is the volume fraction of species A, g is the homogeneous Gibbs free energy of mixing and κ is the gradient energy parameter. The equations are scaled with the following time and length scalings:
where R G is the radius of gyration for the polymer and D AB is the diffusion coefficient. For polymer blends, g is well represented by the Flory-Huggins equation:
and κ is given as follows (Ariyapadi and Nauman, 1990) :
While it is possible to use the default FiPy Viewer to view the simulation results, the use of ParaView enables students to explore more advanced post-processing techniques. The system Gibbs free energy G System as given by eq. (24) was evaluated at multiple time-steps in ParaView and is presented in fig. 6 . It is evident that the Gibbs energy has started to plateau. However, if students continue to let the simulation run, they would be able to observe Ostwald ripening where the smaller blue globules in fig. 6 "disappear" while the larger globules increase in size while G System continues to decrease.
From a numerical perspective, the solution of the modified Cahn-Hilliard equation becomes particularly challenging if students specify excessively large values of N A , N A and/or χ AB . For the symmetric system where N A = N B , typically the simulation diverges at N A χ AB ≥∼ 6. It must be noted that the development of stable schemes for the numerical solution of the Cahn-Hilliard equations is still an active area of research (Chen et al., 2019) .
A total of 9 students signed up for the remote course. We delivered two surveys (referred to as pre-course survey and postcourse survey) to assess the change in student self-efficacy and interest in using Python to solve PDEs. In the post-course survey, we also asked about various aspects of design and delivery of the course. In addition to the surveys, we also held focus group at the end of each session of the course to get further details about the content. The analysis of student feedback and survey data is broken down into three categories: Introduction and code familiarisation, content and structure and lastly course delivery and integration.
The development team took significant care to ensure that the set-up process was as smooth as possible. This involved selecting Anaconda and the Spyder IDE as the environment for students to run code, compiling extensive notes with annotated screenshots for set-up. The student response to comprehensiveness and ease to follow these instructional documents is shown in fig. 7 . While most of the students agreed that the instructions were comprehensive and easy to follow, there were some students who struggled to set up the environment. This was due to older versions of Anaconda which required uninstalling and installing the latest version to avoid errors. To ensure proper set up, the development team held multiple drop-in sessions for those who were facing issues.
Additionally, the first hour of the course was set aside to help student resolve their set-up issues. A common theme in the feedback provided was that more explicit instructions were desired. One student even gave a suggestion of preparing a stepby-step video of the set-up process together with additional information on troubleshooting, indicating it might be more engaging for students compared to written notes.
"If you just had that as like a video and then if people get stuck there can be like Oh OK, so the video tells me I just need to do this so..." (a) Polymer blend morphology @t = 20000. The volume fraction colour bar is included for reference.
(b) Evolution of G System over time, demonstrating how the demixing process decreases the total Gibbs energy. This feedback served as a reminder that even though we had designed the course material to be comprehensive and accessible by our own standards as both educators and users of these research codes, there is always room for improvement. Almost everyone who starts coding will inevitably get frustrated with setting up certain things e.g. an important piece of software or a dependency or even with writing basic code. While such experiences may be instructive, they can be particularly demotivating to beginners. We asked students to rate the activation barriers/difficulties associated with starting and applying Python and the results are shown in fig. 8 . Environment set up, familiarization with syntax, and lac of aim/problem to solve posed high activation barriers to programming according to majority of the students. Therefore, due care is needed to reduce this friction as much as possible to create an enabling environment for experiential learning (Inguva et al., 2018; Shah et al., 2020) .
As most participants did not have any experience in Python, minimal emphasis was placed on syntax with students advised to treat the code as boilerplate, i.e., default codes which they can copy and paste appropriately when needed. We also included the necessary syntax that students would need to be able to adapt the solution as comments that students could comment in or out as desired. This approach worked quite well with most students being able to execute code and make small modifications independently. However, students were generally keen to learn more about Python and how FiPy works and were suggesting that the inclusion of a crash course / refresher in Python might be helpful.
"If it's just copying and pasting, I think it'll be right. I was more thinking in the future. For example, if I want to use this myself to solve some problems. Then it would then be nice to also be introduced to how this whole thing is structured."
Students broadly enjoyed the content of the course, with the visualisation aspect of the numerical solutions being the most well received. This is not a surprising result as providing a visual component can greatly facilitate understanding of concepts, especially more abstract topics such as advanced mathematics (Guzman, 2002) . An exemplar student comment regarding this point is shown below.
"I felt the course explained PDEs better than we got taught it, because you can see it as it was happening and see how one thing affects the other. The graphs gave a more visual approach than when it was taught..."
We were pleasantly surprised when students demonstrated that they were able to keep up with the pace of the course, especially on the second day when more theoretical and complex aspects of numerical methods, PDE theory and thermodynamics related to exercises three and four were introduced. Students broadly acknowledged that the course structure was logical and facilitated their progression through the course material. The further division of exercises two and three into distinct parts with clear goals also helped students by making clear the learning point for that particular part of the exercise.
"I like the fact that you separated into like parts ABC so people could see which code coincided with which sort of change in like the equation that we're making..."
We asked student questions related to their self-efficacy including their self-perceived confidence to a) master and b) apply concepts of PDEs in the pre-course and post-course surveys. The results are shown in fig. 9 . Students who expressed no confidence or low confidence in the pre-course survey showed higher confidence in the post-course survey suggesting the course improved their self-efficacy related to PDEs. The course had most benefit in perception of students who were extremely unlikely to use programming to solve PDEs. One of the key objectives of the course was to provide a starting point for students to use Python for solving PDEs and several students responded favorably to using the codes as templates in the future ( fig. 9 ).
Due to the current COVID-19 pandemic, the course had to be run remotely. Fortunately, due to the various features made available in online collaboration platforms such as Microsoft Teams, the course could be delivered smoothly with no students indicating that remote delivery adversely impacted their experience. One way Microsoft Teams facilitated this process was through the feature of being able to create polls in the chat function. This enabled teaching staff to create polls at specific intervals to gauge students' understanding of a specific point or to gather information on what aspects of the code did students independently explore during breaks.
Most students indicated a desire for more active and structured learning components related to the actual code and the running of simulations during the course. The original design of the course centered around providing students the flexibility to independently explore various aspects of each exercise such as different boundary and initial conditions or different numerical schemes with the exercise script serving as a basis for students to then modify. This approach was not favoured by students. Many of the suggestions regarding this point indicated that students desired a couple of formal and structured questions or tasks which they could execute at the end of each exercise and compare to a correct solution to consolidate their learning. We intend to implement this as part of the second iteration of the course.
"If you had a couple of exercise questions or stretch questions for your knowledge at the end and you were like we would like you to obtain this graph or we would like you to change this parameter and think about the understanding that it has or model this question(from the problem sheets or exams) ... Something like that, and then you can have a predictable outcome that people can reproduce with just a few exercise questions at the bottom I think would be helpful."
Regarding the integration of the activities from the course into lecture based courses such as the Math II module, the first and most readily adaptable aspect is the generation of high quality plots and animations for visualisation purposes. Beyond that, the inherently open-ended nature of such a tool lends itself easily to being used as the basis of a coursework which could be based on a specific example from the range of equation systems covered in the present work or even a new system that could be of interest to students.
We have demonstrated how the PDE solver code FiPy can be effectively introduced to students, with very little coding experience, in a short period of time. The key to this is designing the course to be as accessible as possible and design-(a) Student response to their self-perceived confidence to master the concepts of PDE.
(b) Student response to their self-perceived confidence to apply the concepts of PDE.
(c) Student response to their likelihood of using programming to solve PDEs (pre course vs. post course).
(d) Student responses to their confidence in solving problems and using Python in the future after the course. ing each exercise to build on previously taught material, enabling faster progression. Referencing the four factors outlined previously for evaluating codes for use in the classroom, FiPy and other PDE solver codes in Python do well in all four aspects. The open-source nature of all the tools used in the present work: Anaconda, FiPy and ParaView, mean that availability is not an issue. These tools are also frequently used in industry and academia, demonstrating their scalability. Lastly, we have demonstrated how a variety of problems of different complexities can be easily implemented and rolled-out to students, highlighting the ease-of-use of these tools and the ability for them to be tailored to meet a variety of learning objectives. This course improved student perception to Python and their self-efficacy in concepts of PDEs with students with less confidence showing the greatest boost.
The present work also reveals various multiple challenges and issues educators should take note of when rolling out such an introductory course. Despite our best efforts to make the set-up process as seamless as possible, students indicated this still remained a significant pain point. Fortunately, receiving such student input at a small scale course means that this can be addressed in future iterations. Regarding the structure of the course, educators should provide more structured tasks during the content delivery process so that students can consolidate their learning and remain engaged. These tasks could take the form of simple tasks such as modifying the initial/ boundary conditions of a problem, changing values of certain dimensionless parameters, and comparing it to a model solution.
The highly flexible nature of these codes enables them to be used in a variety of settings beyond the introductory course we have presented. A review of the various user discussion forums or code documentation reveals several more capabilities and use cases that can be integrated into the classroom in a variety of topics at different academic levels. We hope to explore the nuances between the different codes and discretisation schemes and its integration into a more advanced course of numerical methods in future work.
Authors would like to thank students who participated in this course and provided their valuable feedback.